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Foreword 


Safran Tech is the corporate research center of Safran Group, a French multinational 
company that designs, develops and manufactures aircraft and rocket engines, as 
well as various aerospace and defense-related equipment. This book is the result of 
the collaboration of two of my research engineers at SafranTech, Fabien Casenave 
and Nissrine Akkari, with David Ryckelynck, professor of Applied Mathematics at 
MINES Paris PSL, a French engineering school. This collaboration dates back to the 
creation of Safran Tech in 2015 and has involved Ph.D. students, namely Thomas 
Daniel and Hamza Boukraichi. Thomas Daniel won the 2022 AMIES Ph.D. award 
(AMIES: Agency for Mathematics in collaboration with companies and the society) 
for his work on dictionary-based ROM-nets, which are introduced in this book. 

Nowadays, the engineering and design processes involved when creating a new 
mechanical part, component or product in industrial companies like Safran rely on 
some level of numerical simulation. In some cases, complex and high-dimensional 
models are required to be computed, and often in a so-called “many queries” context: 
when searching for an optimal design, quantifying uncertainties or exploiting a “real- 
time” simulator, the physics problem is solved a large number of time at different 
configurations. Hence, we need fast surrogates, with a control of the approximation 
error. Among the available families of methods, this book deals with projection- 
based reduced-order models, which consists in solving the same physics equations 
as the reference high-dimensional models, but on a reduced dimensional vector 
space. Resorting to the known physics equations provides explainability, error esti- 
mates in some cases, and good performance with sparse learning databases. Clas- 
sical and modern (e.g., deep learning) machine learning technologies are used to 
assist projection-based reduced-order models, either in annex tasks like clustering 
or classification or by providing non-linear dimensionality reduction procedures. 
Challenges of these methods include geometrical variabilities, predictive uncertainty 
quantification, and efficiency (i.e., speed-up factor, also taking the learning stage into 
account). 

This book is addressed to graduate students and researchers who look for an 
introduction to projection-based reduced-order models, an understanding of the main 
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features and challenges of these methods, as well as its state of the art application to 
real-life engineering design problems and its extensions. 


Paris, France Christian Rey 
September 2023 
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Chapter 1 A) 
Structured Data and Knowledge in geai 
Model-Based Engineering 


1.1 Nomenclature 


In the book, the following notation is used: bold face symbols denote vectors (low- 
ercase letters) a, b ... or matrices (uppercase letters) A, B... Data are denoted by 
X for input data and labels or output data are denoted by Y. The spatial gradient 
operator is denoted by V, and the divergence is expressed by V - e. The dependence 
on spatial coordinates & is most of the time omitted for simplicity of notation. Mod- 
eling parameters are denoted by u or pg when they have a tensor shape (multiple 
indices). Weights in neural networks are denoted by W. In mechanical models, the 
Cauchy stress is denoted by o.. Solutions of physics-based models are denoted by ù 
for partial differential equations (PDEs) or ordinary differential equations (ODEs). 
High-fidelity solutions in finite dimensional space are denoted by u. Solutions in 
reduced dimensional spaces are denoted by t. PDEs are set up over a domain denoted 
by 2. Solution spaces are the following: 


e V for PDEs, 
e V, for high fidelity solution space of dimension N, 
e V, for reduced vector space of dimension n < N. 


The finite element shape functions are denoted by ø. The coordinates of a finite 
element solution are denoted by the vector u. The reduced vector space V, is spanned 
by empirical modes denoted by (W%)x=1,...n. The matrix form of a reduced basis is 
denoted by V € R%*". The vector of reduced coordinates is denoted by y € R”; 
Yk, k € [1,...,n] is the kth component of y. When needed, variables depending on 
parameter yw are denoted using an exponent è” or an index e,. 
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1.2 Model-Based Engineering 


Today, a large set of engineering tasks is supported by mathematical models related 
to various scientific disciplines. This set of tasks is called Model-Based Engineer- 
ing (MBE). In this book, we restrict our attention to geometrical or morphological 
models, thermal models, mechanical models and statistical models, at various scales. 
The related mathematical equations can be algebraic equations, ordinary differential 
equations, partial differential equations (PDEs), or every possible combination of 
these types of equations in a coupled system of equations. 

Various numerical methods have been developed in applied mathematics, to find 
approximate solutions of these different types of equations. Therefore, mathematical 
properties are available related to the convergence of the numerical solutions to the 
exact solutions having no numerical approximation. The engineers have a strong 
confidence in the numerical predictions obtained by these models, provided that 
they are well used in their domain of validity. Hence, in practice, numerical models 
and numerical predictions are ubiquitous in MBE. 

In industrial sectors, such as automotive industry, aeronautical industry, energy, 
shipbuilding, electronics manufacturing, etc., computer platforms for MBE have 
been developed for assemblies of components, to ease the communication of data 
and models, between engineering tasks. Quality and efficiency in engineering tasks 
depend directly on the continuity of data on computer platforms. As explained in 
the ATLAS program of AFNet,! such continuity of data requires the development of 
standards for data exchanges between stakeholders in a project, a production line, a 
department, etc. Assemblies of components that are equipped with numerical models 
are termed complex systems in this book. 

A complex system is an assembly of components that have a numerical represen- 
tation related to coupled physics-based models, for its design phase, its manufacturing 
phase or its exploitation phase. 

In computer platforms for MBE, complex systems have an idealized digital rep- 
resentation, termed digital mock-up, where model couplings can be implemented 
in a numerical model dedicated to an engineering task. For instance, in continuum 
mechanics, the mechanical model of a component is weakly coupled to the geomet- 
rical model of this component. Many multi-physics couplings in coupled problems, 
can be implemented on computer platforms for MBE. Complex numerical simu- 
lations in MBE are supported by computers that range form laptops to computer 
clusters for high-performance computing (HPC). In most cases, geometrical models 
are created by using computer-aided design (CAD) in order to obtain a common 
geometric model for both design tasks and manufacturing tasks, of each part and 
assemblies in a complex system. 

Data flux in computer platforms for MBE are related to design tasks, manufac- 
turing tasks, non destructive testing (NDT) tasks, exploitation or monitoring tasks, 
inspection tasks, repairing tasks, recycling tasks. Figure 1.1 is a summary of these 
data fluxes. 


l https://www.afnet.fr/en/. 
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Fig. 1.1 Schematic view of data continuity in an computer platform for Model-Based Engineering 
(MBE), enabling engineering tasks related to CAD, PDE solution 
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The value of data, beyond the tasks at hand, should not be ignored by engineers. 
In particular, the value of data must be considered in light of advances in artificial 
intelligence (AI) and machine learning methods. 


1.3 Motivation 


Our main paradigm is: the more we know about underlying mathematical properties 
or physical properties of data, the less we need data for machine learning tasks. In this 
book, we aim to learn manifold from data and conventional knowledge in engineering. 
This paradigm has been recently employed in Physics Informed Neural Networks 
[6, 7] and knowledge transfer [4]. Similarly to transfer learning, we are “motivated 
by the fact that people can intelligently apply knowledge learned previously to solve 
new problems faster or with better solutions” [5]. 

Much data in computer platforms for MBE are structured data, in the sense that 
their format is predefined and documented (this notion is different from structured 
and unstructured meshes, referring respectively to constant rectilinear and hetero- 
geneous meshes). These structured data are supported by a graph or they have a 
tensor format, like digital images. Since deep learning has made huge progresses in 
computer vision, image transformation and for graph neural networks, there is an 
opportunity to develop deep learning in MBE. The development of sensors and image 
acquisitions in non destructive testing of complex systems makes this opportunity 
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more and more relevant. The connection of computer platforms to HPC facilities 
leads to the two following issues: 


e Where are the data in MBE? 
e What can be learned that is useful or valuable to engineers? 


In this book, we propose numerical methods and algorithms to assist engineers in 
various conventional tasks by using machine learning. This gives value to data that 
flow in computer platforms for MBE. Our goal is to augment engineers capabilities 
by using artificial intelligence, without ignoring the knowledge already available 
for engineering in all scientific disciplines. We think that the collaboration between 
engineers and data scientists is going to grow up, but engineers will not be replaced 
by them. At the engineering school MINES Paris PSL, professors David Ryckelynck 
and Elie Hachem created in 2017 master courses to ease such collaborations between 
engineers and data scientists. The title of this set of courses is Ingénierie Digitale des 
Systémes complexes (Digital engineering of complex systems). The creation of such 
a course would not have been possible without the help of SAFRAN, both for the 
illustrations of classical AI and for the development of an engineering augmented 
by a hybrid AI, which preserves scientific knowledge. 

Acceleration of engineering tasks is the main objective here, but we are also able 
to develop engineering tasks that were unaffordable without machine learning. This 
is especially true when considering uncertainty quantification (UQ) or image-based 
digital twining of complex systems. 

Following the definition in [2], a digital twin involves similar numerical models 
to mock-up but it is related to an existing instance of a complex system or an instance 
of one of its component. A digital twin also include the related observational data. 
Hence, in essence, a digital twin has multimodal data, such as observational data 
and simulated data. The acceleration of image-based modeling for digital twining of 
manufactured parts is going to foster the development of functional assessment, in 
which defects harmfulness are evaluated with respect to physics-based of mechanical- 
based specifications. 

The knowledge accumulated via digital twining or via sensors only, enable the 
development of uncertainty quantification. Sensors and numerical mock-up being 
more and more spatially resolved and accurate, the acceleration of numerical predic- 
tions is becoming crucial for practical application in engineering. In Sect. 5, mechan- 
ical predictions are accelerated by using machine learning in order to account for 
stochastic variation of temperature fields when computing the lifetime of turbine 
blades. 

In this book we develop the following machine learning tasks, via manifold learn- 
ing: 


classification, 

anomaly detection, 
regression, 
dimensionality reduction. 
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Data visualization is out of the scope of this book, although dedicated manifold 
learning techniques have been proposed in the literature. In this book, manifold 
learning is thought of as an attempt to generalize linear frameworks for numerical 
methods in engineering. 

When accumulating data in computer platform for MBE, data do not fill the entire 
ambient space used for the encoding of these data. There exists an empirical space, 
or a latent space, occupied by data whose dimension is usually much smaller than 
the ambient space. Dimensionality reduction and manifold learning aim to find a 
numerical approximation of this latent space. As models are ubiquitous, we have 
to consider large sets of synthetic data, or more precisely, sets of simulated data 
for engineering tasks. These simulated data are solutions of numerical equations. 
They belong to a solution space, which is the ambient space for simulated data. A 
solution space is directly related to the general numerical scheme used to compute 
instances of numerical predictions. Given an engineering task, an ambient space is a 
common space for all data instances. It does not depend on any variable parameter. 
The dimension of the ambient space is set up so that each instance of data is a point 
inside the ambient space. 

Two kinds of latent spaces can be defined when considering simulated data: 


e latent space for surrogate modeling, 
e latent space for projection-based model order reduction. 


Projection-based model order reduction is detailed in Chap. 2. Our recommendation 
for selecting a surrogate model rather than a projection-based reduced order model 
is the following. If you consider a mature engineering task so that: 


e you can decide a priori quantities of interest, 

e all influencing parameters have known range of variation, and if possible a known 
probability density functions, 

e you are able to sample the parameter space without suffering from the curse of 
dimensionality and build a data-base, 

e you need real-time predictions, 


then you should use a surrogate modeling rather a projection-based reduced order 
modeling. But, when modeling a complex system by using PDEs, the selection of 
simulation outcomes is rarely trivial. Generally speaking, projection-based reduced 
order models are slower than surrogate models, but one can adapt the outcomes 
of reduced predictions and their validity domain is not restricted to the parameter 
domain used to learn the reduced order model, both in term of extent and in term of 
dimension. 

In practice, data accumulation has to be supplemented by data compression 
schemes but also data pruning scheme. As explained in [3], manifold learning for 
projection-based model order reduction is a convenient way for data pruning that 
preserves modeling capabilities. Pruning algorithms facilitate data augmentation for 
high-dimensional data, by limiting the memory requirements for augmented data [1]. 
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So, data in MBE are simulated data and observational data, where simulated 
data are synthetic data forecast by physics-based models. Unfortunately, observa- 
tional data on complex systems are extremely rare in engineering. For instance, air- 
craft landing gear tests, automotive crash tests, etc. are mostly done for conformity 
assessment processes or for certification processes. For example, in the aeronautics 
industry, aircraft manufacturers must demonstrate compliance of new products with 
regulatory requirements. This compliance demonstration is done by analysis during 
ground testing (such as tests on the structure to withstand bird strikes, fatigue tests 
and tests in simulators) but also by means of tests during flight. Usually, manufac- 
turers have many more instances of observational data, at the components scale, 
during manufacturing processes, maintenance and exploitation tasks. In engineer- 
ing, the product lifecycle management (PLM) refers to the management of data and 
processes across the entire lifecycle of a component or a product. PLM requires a 
high level of data continuity across the supply chain, by developing data standards 
and dedicated computer platforms. Such computer platforms open the route for more 
machine learning in engineering, especially in MBE where engineers are facing very 
complex structured and unstructured data. 

In this book we restrict our attention to spatially structured data. The detailed 
description of the data we worked with are available in Chap. 6. 
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Chapter 2 A) 
Learning Projection-Based cieca; 
Reduced-Order Models 


2.1 Motivation and Basic Assumptions 


In this chapter, we introduce the solution space for high-fidelity models based on 
partial differential equations and the finite element model. The manifold learn- 
ing approach to model order reduction requires simulated data. Hence, learning 
projection-based reduced order models (ROM) has two steps: (i) an offline step for 
the computation of simulated data and for consecutive machine learning tasks, (ii) an 
online step where the reduced order model is used as a surrogate for the high fidelity 
model. The offline step generates a train set and a validation set of simulated data. 
The accuracy and the generalisation of the reduced order model is evaluated in the 
online step by using a test set of data forecast by the high-fidelity model. The test set 
aims also to check the computational speedups of the reduced-order model compare 
to the high-fidelity model. 

Learning projection-based reduced order model makes sense only if there is a 
significant computational speedup at the price of an acceptable loss of accuracy in 
predictions. The longer the computational time of the high-fidelity model, the smaller 
the acceptable speed up, if we save hours or days of numerical simulations. Regarding 
the acceptable accuracy of reduced predictions, engineers working in industry and 
scientists working in laboratories do not have the same expectations. In our own expe- 
rience, learning projection-based reduced order model has the capability to adapt to 
engineering tasks, high-fidelity models elaborated in laboratories, in terms of accu- 
racy and computational time. This approach contributes to data continuity between 
physics-based knowledge developed in laboratories and practical applications for 
engineering tasks. 

From the mathematical point of view, Céa’s lemma gives an overview of manifold 
learning for model order reduction applied to elliptic equations. Let V be a real 
Hilbert space, such that the weak form of the elliptic equation reads: find X € V, 
such that: 

aŭ, d) = LÆ, VUEY, (2.1) 
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where a(., -) is a bilinear form, with coercivity constant 6 > O and continuity conti- 
nuity constant C“ > 0. L(-) is a linear form. Section 2.2 gives more details on weak 
forms. Let us V, a finite dimensional subspace of V. Here, V, is an approximate 
solution space for the elliptic equation. The approximate solution of the elliptic equa- 
tion is denoted by u € Van C V. The Galerkin projection of the elliptic equation onto 
the solution space reads: find u € V} C V such that: 


a(u,v)=L(v), Vue Vy, (2.2) 


where V, has been substituted for V. Céa’s lemma states that: 


a 


a CF ik 
lu —ul| < T lu — vl, Vues, (2.3) 


where || - || is a norm in V. The related scalar product is denoted by (-). In Finite 
Element models, V, is the span of finite element shape functions. But Céa’s lemma 
holds for all finite-dimensional subspace of V. The closer % to the solution space V}, 
the smaller the right-hand term of Céa’s lemma (2.3), and therefore one can expect 
a better prediction u in Eq. (2.2), although we do not know %. 

In few words, a reduced-order model is obtained by introducing a smaller solution 
space V, C Vp of smaller dimension n < N. A projection-based reduced order 
model can be achieved by using the Galerkin projection (2.2), where V, is substituted 
for V;,. The conclusion of Céa’s lemma holds again. The closer Ù to the solution 
space V,, the better the prediction # € ‘V,,. Manifold learning comes into play when 
we are given a set of predictions (u);—1..m in a common ambient space V, related 
to a given finite element mesh. The basic assumptions in manifold learning for model 
order reduction are: 


e a latent space of reduced dimension is hidden in the data (u®);=1 
is denoted by n, 

e a machine learning algorithm is available to learn this latent space by using a train 
set of simulated data extracted from (u®);=1.....m, 

e the distance between u and this latent space is small enough, although we do not 
know iw, 

e a numerical scheme enables the projection of the elliptic equation onto the latent 
space, in order to set up the reduced order model, 

e the computational complexity of the solution of the reduced order model is smaller 
that the computational complexity of the finite element prediction, 

e the computational complexity of the reduced order model is an increasing function 
ofn. 


m» its dimension 


As explained above, when the latent space is a vector subspace V,, both Galerkin 
projection and Céa’s lemma hold, but simulation speedup may not be achieved. The 
study of more complex situations is the purpose of this Chapter. An estimation of 
computational complexity of projection-based reduced order model is proposed in 
Sect. 2.3.6. 
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Remarks: 


In the Rayleigh-Ritz method a small set of trial functions that satisfy the boundary 
conditions for % is introduce to span a solution space. This solution space is not 
related to any finite element model. The inclusion of the latent space in the ambient 
space V, is essential for model order reduction. 

The finite element ambient space V, incorporates homogeneous Dirichlet bound- 
ary conditions (Ù = 0) on a boundary of the domain where the partial differential 
equations are set up). When such boundary conditions may change in the instances 
aË) i=1,...,m these conditions must be taken into account as a linear constraint that 
supplement the partial differential equation. Such an issue appears when consid- 
ering contact problems [38] for instance. 

Important limitations of projection-based model reduction methods includes sit- 
uations where the geometry has to be handled in the exploitation phase of the 
reduced-order models, for instance when the problem features contact boundary 
conditions, crack propagation or when the geometry is a variability of the problem 
to learn. Geometrical variabilities are handled in the authors’ works [1, 2, 22, 60, 
61, 92]. 


2.2 High-Fidelity Model (HFM) 


Consider an abstract partial differential equation in a domain Q, with a y-variability: 
D, &,p)=0, EEQ, WEY. (2.4) 


The weak form of this partial differential equation reads: find 7 € V such that 
[ wom. wag =o. Woe. (2.5) 
2 


As an illustration, the concepts of this chapter are illustrated on a nonlinear struc- 
tural mechanics problem, for which details on the high-fidelity model are provided in 
this Section. For an example with another physics, we refer to the authors’ work [21], 
where a nonlinear transient thermal problem is considered. 

The mechanical structure occupies the domain Q, whose boundary 0{2 is parti- 
tioned as 0Q = Np U OQ such that IQ N IQS, = Y, see Fig. 2.1. 

The structure is subjected to a quasi-static time-dependent loading, composed of 
an homogeneous Dirichlet boundary conditions on 0Q2p and Neumann boundary 
conditions on dQ, in the form of a prescribed traction Ty, as well as a volumic 
force f. The setting depends on some variability u, which can be a parameter vector, 
or represent some nonparametrized variability. The evolution of the displacement 
UF (E, t) over (£, t) € Q x [0, T] is the solution of the following equations: 
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Fig. 2.1 Schematic dQNn 


representation of the è 
structure of interest [18] r we B 


a 
öp 

Ey 1 = a : paet ; 
eQ) = 5 (var + (varT) , in Q x [0, T], (compatibility equation) 
div Gi ) +f" =0, in Q x [0, T], (equilibrium equation) 
oÉ =o-(e(u"), y”), in Q x [0, T], (constitutive law) 
a” = 0, in dQp x [0, T], (prescribed zero displacement) 
ob -NgQ = TE, in Qy x [0, T], (prescribed traction) 
qr = 0, y# = 0, in Qatt = 0, (initial condition) 


(2.6) 
where € is the linear strain tensor, o! is the Cauchy stress tensor, y“ denotes the 
internal variables of the constitutive law and nag is the outward normal vector on 
dQ. We precise that the evolution of the internal variables y” is updated when the 
constitutive law is solved. 

Define Hj (Q) = {v € L7(Q)| n € L?(Q), 1 <i <3 and vlag = 0}. Denote 
{@ihi<i<n € R** | a finite element basis whose span, denoted V,, constitutes an 
approximation of Hj (Q)?; N is the number of finite element basis functions, hence 
the number of degrees of freedom of the discretized prediction. A discretized weak 
formulation reads: find u” € V, such that for all v € Vp, 


f o(e(u"), y) : €(v) = i Jeu +f TU; (2.7) 
Q Q IQ 


that we denote for concision F” (u“) = 0, where u” is the vector of N coordinates for 
u” € V. A Newton algorithm can be used to solve this nonlinear global equilibrium 
problem at each time step: 


oe (ut) (at — ult) = FH (uh), (2.8) 
where 
DF" |, (eet 
D (u') = f ele) (ew), 9") : €), (2.9) 
u ij 2 


and 
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Linear data compression 


u u(y) 
Variability with small influence on u « Small » solution manifold 
Parameter space Solution space # c L? (Q) 


Fig. 2.2 Linear manifold learning 


F“ = fo (eut, y”) TOR EA fe Ty pi. (2.10) 


In the two relations above, K (e (u5), yë ) is the local tangent operator, u“" € V 


is the k-th iteration of the discretized displacement field at the current time-step, and 
N 


uk = (u) € RY is such that u”* = > u!" pi. In particular, f”, TH, wt* 
I<i<N = 

and y“ are known and enforce the time-dependence of the solution. Depending on the 

constitutive law, the computation of the functions (ah, y) > e (e (ul), y“) and 

(u, y") K (e (ul*), y“) can require the resolution of a complex differential- 

algebraic system of equations. 


2.3 Linear Manifold Learning for Projection-Based 
Reduced-Order Modeling 


We start by explaining the online phase. Since we want to construct and solve the 
reduced-order model (ROM) in the most efficient way, the offline phase is dedicated 
to precompute as many steps as possible, under the considered variability. 

Linear manifold learning means that the solution manifold is approximated by a 
vector subspace of the ambient solution space, as illustrated in Fig. 2.2. 


2.3.1 Approaches Preceding the Use of Machine Learning 


In structural mechanics, normal modes have been introduced for the analysis of 
vibrations in structures. When considering free vibrations, without external force, the 
solution of linear hyperbolic equation is sought by using the separation of space and 


14 2 Learning Projection-Based Reduced-Order Models 


time variables: uy (x, t) = w(x) u(t), where x € Q is the space variable and t € R 
the time variable. The hyperbolic equation of free vibration reads: find w(x) € Vy 
and u(t) € R such that uy (x, t) = Y (x) T(t) and 


(p tin, Vy) ta(uy, Yy) =0, Vvy € Vy. (2.11) 


It follows that T(t) is an harmonic function of frequency denoted by f and w is 
the eigenvector related to the eigenvalue À = (2 2 f)? of the following generalized 
eigenproblem: find y € Vy and A such that 


a(W,vv) —A (pW, vw) =0, Vw € Vy, (2.12) 


where the rank of this system of equation is supposed to be N — 1 in order to find non 
zero eigenmodes. This eigenproblem admits N orthogonal normal modes that span 
Vy. Therefore a selection of n normal modes span a reduced subspace of dimension 
n. In the beginning of the 21st century, model reduction using variable separation 
scheme in partial differential equations has been extended to an arbitrary number of 
variable by using low-rank approximations such as the Proper Generalized Decom- 
position [5]. Eigenmodes are global functions in contrast to finite element shape 
functions that have a local support. For dynamical problems involving nonlinear 
contributions to the PDE, adaptive computations of reduced subspaces have been 
proposed by Almroth et al. [4] and Noor et al. [79], by using Rayleigh-Ritz global 
functions. The set of these global functions is a reduced basis of the finite-element 
approximation-space. 

The idea of using statistics to generate a solution space for differential equations 
has been proposed in the seminal work of Lorenz in [69] (in page 31), by using 
empirical orthogonal functions. The Galerkin projection of PDEs on empirical modes 
have been first developed in [13], where a reduced basis is computed via the proper 
orthogonal decomposition [70] of observational data. This was the first step towards 
manifold learning for projection-based model order reduction, which we now present. 
For other presentations of this technologies, the reader can refer to [81, 87]. 


2.3.2 Online Phase: Galerkin Projection 


The reduced-order model is constructed in the form of a Galerkin method written 
on a Reduced-Order Basis (ROB). In the present case, it consists in assembling the 
physical problem in the same fashion as the HFM in Sect. 2.2, with the difference that 
the finite element basis (¢;) |<< E€ RY*™, is replaced by a ROB (Wj) 1<j<n E€ R"™™, 
with n < N. Hence, the reduced Newton algorithm is constructed as 


DF 
Di 


(a) (ye! = yi) =—-Fy (ii',) i (2.13) 
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where 


DF un /. 
A (â) 


= (wi) : K (€(a@4), yu)  € Wa) (2.14) 


E 
ij Q 


and 
F, (at), = f o (eai) sew- f fhn- f Tae (215) 
Q Q dQn 


In the two relations above, âk cV := Span (W;)i<i<n is the k-th iteration of the 
reduced displacement field for the current time-step and y“* = (vf ) e R” 
l<i<n 


is such that 


it = > y Wy. (2.16) 
i=l 


Notice that the use of the Galerkin method is made possible by the linearity of the 
tangent problem (2.13) and the choice of a linear dimensionality reduction technique 
in (2.16). 

The online stage is called efficient if the reduced problem can be constructed 
and solved in computational complexity independent of N. When the variability u 
is parametrized, efficiency is possible by precomputing various terms. With non- 
parametrized variability, depending on its nature, some assembling task with a linear 
complexity in N may be required at the beginning of the online stage (for instance 
for a boundary condition). All these scenarios are handled by genericROM, the ROM 
library developped at Safran and presented in Sect. 4.2. 

The offline phase contains three steps, for which we present below the method- 
ological choices made in genericROM. 


2.3.3 Offline Phase 
2.3.3.1 Data Generation 


This step corresponds to the generation of the snapshots by solving the high-fidelity 
model. In parametric contexts, the simplest workflows consist in choosing param- 
eter values a priori, following Design of Experiments (DoE, see [50] for a recent 
technique), and computing the corresponding snapshots by solving the high-fidelity 
model. 
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2.3.3.2 Data Compression: Dimensionality Reduction 


This step corresponds to the generation of the ROB (W;)1<;<,. One of the most classi- 
cal and simple method is the snapshot Proper Orthogonal Decomposition (POD) [25, 
86], detailed below: 


1. Choose a tolerance €pop. 

2. Compute the correlation matrix C; j = Ís ui- uj, l <i, j < Ns, where N, is the 
total number of HFM snapshots. 

3. Compute the epop-truncated eigendecomposition of C: & € R% and A; > 0, 


where 1 < i < n, are the n first eigenvector and eigenvalues. 
N, 


1 s 
uj(x)é&;,1<i <n. 

The advantages of the snapshot-POD are a reasonable computational complexity 
when the number of degree of freedom of the high-fidelity model are much larger than 
the number of snapshots, and the fact that this algorithm can be easily parallelized. 

Variants can be used, for instance in the Reduced Basis [84] and the POD- 
greedy [41] methods, with respectively an orthonormalization of the computed high- 
fidelity snapshots and the incremental Singular Value Decomposition (SVD) [15]. 


4. Compute the reduced order basis W(x) = 


2.3.3.3 Operator Compression 


A ROM is called online-efficient if in the online stage, the reduced problems can be 
constructed and solved in computational complexity independent of N. The operator 
compression step consists in additional treatments required for the efficiency of the 
online stage, by pre-processing some computationally demanding integration tasks 
over the high-fidelity domain Q and dQ y. We notice that without any additional 
treatment, the numerical integration involved in the assembling of Eq. (2.13) strongly 
limits in practice the efficiency of the ROM: no speedup with respect to the high- 
fidelity model can be obtained in practice. The complexity of such an additional 
treatment depends on the type of parameter-dependence of the problem. This step is 
actually needed for all classes of problems reduced by projection-based methods. 

Consider the simplest case: a linear problem with an affine dependence in the 
parameter u, for instance A „u = c, where A, = Ag + “Ay. Denote V, the matrix 
whose columns are the vectors of the ROB evaluated at the high-fidelity degrees 
of freedom. The obtained ROM writes V” A, Va = V” c: it is not assembled in the 
online phase, but rather the matrices VT AgV and VT A,V and the vector V"¢ are 
precomputed in the offline stage so that the reduced problem is constructed without 
approximation and efficiently by summing two small matrices. The operator com- 
pression step consists in this case in the construction of V” AoV and V7 A,V in the 
offline stage. 

Actually, there exist linear problems for which the operation compression step 
require an additional approximation, and nonlinear problems that can be carried- 
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out exactly. In the first case, consider A,,u = c with Ajj = 1 V (g(x, uj (x)) . 
Q 


Vo;(x) and c; = Íe f(x)pi(x), where u is the unknown, f a known loading and 
g(x, u) a known function with variables x and u that cannot be separated: the 
previous precomputation of reduced matrices cannot be applied, and a treatment is 
required to, for example, approximately separate the dependencies in x and u of g 
as g(x, m) = Y_i ge(x)ge(w). Then, VA V ~ Yi, gt (u)Ar where (Ax);; = 


f V (sf (x); (x)) - Vg; (x), so that the efficiency of the online stage is recovered; the 


Empirical Interpolation Method has been proposed in [14, 71] for this purpose. A case 
of linear problem in harmonic aeroacoustics with nonaffine dependence with respect 
to the frequency is available in [20]. Conversely, nonlinearities can be handled without 
approximation in some cases, for instance, the advection term in fluid dynamics can 


be precomputed in the form of an order-3 tensor: | Wi- (Wr; . V) Wiel <i, j,k <n; 
Q 


see [3] for the reduction of the nonlinear Navier-Stokes equations with an exact 
operator compression step. Other examples are found in structural dynamics with 
geometric nonlinearities, where order-2 and -3 tensors can also be precomputed, 
see [58, Sect. 3.2] and [73]. 

When additional approximations are required, the methods proposed for the oper- 
ator compression step are call “hyper-reduction” in the literature. This term was 
coined by the seminal method proposed in [88] in 2005, but has been extended to 
refer to all the methods proposing a such second reduction stage. Hyper-reduction 
methods include the Empirical Interpolation Method (EIM, [14]), the Missing Point 
Estimation (MPE, [12]), the Best Point Interpolation Method (BPIM, [78]), the 
Discrete Empirical Interpolation Method (DEIM, [23]), the Gauss-Newton with 
Approximated Tensors (GNAT, [17]), the Energy-Conserving Sampling and Weight- 
ing (ECSW, [36]), the Empirical Cubature Method (ECM, [45]), and the Linear 
Program Empirical Quadrature Procedure (LPEQP, [100]). The reader can find an 
algorithmic comparison of the Hyper-Reduction and the Discrete Empirical Interpo- 
lation Method for a nonlinear thermal problem in [39]. A particular focus is given 
on hyper-reduction techniques via oblique projection and empirical cubature in the 
following sections. 


2.3.4 Hyper-Reduction via a Reduced Integration Domain 


Hyper-reduction via a reduced integration domain has been proposed in [88]. It 
requires a train set of displacement predictions, so that a reduced approximation 
vector space can be trained. The finite elements simulations that generate the train 
set of displacement fields are also predicting stresses fields ø and internal vari- 
ables y. Hence, additional reduced bases can be trained for these variables, by using 
these simulation results [89]. Heuristically, we found it more accurate to include a 
reduced basis for the stresses ø , as an additional reduced basis for this hyper-reduction 
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scheme. Such a reduced basis is also very convenient for error estimation [91]. The 
finite element shape functions for displacement fields are denoted by (@;)j=1,....x. 
For stress fields, we also need to introduce a related finite element representation. 
We denote by (7 );=1,...,ve the dedicated shape functions. In the linear framework 
of manifold learning, we assume that the same finite element mesh is used for the 
target simulation, in the online step or for the test set of data, and all simulations 
used to generate the train set of data. 

In practice, the implementation of the hyper-reduction follows the manifold learn- 
ing step that trains reduced bases for displacements and stresses. We recall that they 
are respectively denoted by V € R‘*” and V° e RM°*"? in their matrix form, and 
(Wy )kat,...n (WY )k=1.,....n7 in their continuous form: 


N 

VO =) 9,@Vie, VxEQ, k=1,...,n, (2.17) 
i=1 
ee 

Ww =) OO VEER, b= Laan’. (2.18) 
i=1 


The reduced displacement reads: 


n 


f(x) = 5 PLX) yk, VX E Q. (2.19) 


k=1 


The hyper-reduction method proposed in [88] aims at computing reduced coordinates 
(Yq) k=1,....n introduced in Eq. (2.19), by projecting the equilibrium equation on V, via 
a restriction of the domain Q to a Reduced Integration Domain (RID) denoted by 
Qr. By following the empirical interpolation method [14], interpolation points are 
computed for column vectors in V and V°” separately [48]. The set of respective 
interpolation points are denoted by P” and P° . We follow Algorithm 2.1, proposed 
for the Discrete Empirical Interpolation Method [23]. 


The RID Qp is such that it contains the interpolation points related to V and V°. 
For engineering applications, the RID can also include a zone of interest in Q that 
is user-defined, by using a subset of finite elements. This zone of interest is denoted 
by Qz; C Q. By construction, for contact-free problems, the RID is the following: 


Qr = Qz Ujepu Supp(Q;) Uiep> supp(Q? ), (2.20) 


where supp( f) € Q is the support of function f. In practice, Qg has its own finite- 
element mesh. It is a reduced mesh involving much less elements than the original 
finite element mesh of Q. It can be enlarge by adding a layer of connected element in 
the original mesh. Integration of constitutive equations are performed on this reduced 
mesh, without any intrusive operation on the original finite element solver. A similar 
hyper-reduction scheme has been developed in [38] for contact problems. 
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Algorithm 2.1: Interpolation points of the Discrete Empirical Interpolation 
Method (DEIM) [23] 
Input : reduced basis vectors V[:, k] € RY k=1,...M 


Output: interpolation point index set J := {i1,...,im} 
1 set Jp := 0; // initialization 
2 for/=1...,Mdo 
3 U <— Vi. 1: d-1)]; // truncated basis 
4 | A} Or Ue, DU I" ; // projector 
5 q: < VI:, 1] — U1 AV[7-1, L] ; // interpolation residual 
6 iy <— arg maxie{1,...,N} (qu)il ; // maximum of residual 
7 h :=h-1 Ui}; // extend interpolation points 
8 end 
9 set I := Im. 


Once the RID is obtained, a set of test functions is set up in order to restrain the 
balance equations to Q. They are denoted by Y% p ;: 


P= fiet. Nh @d2=0}, (2.21) 


Q\QR 
Wri) = >> 9, (x) Vii, Wee Q, j=1,...,n, (2.22) 
icP 


where P is the set of all degrees of freedom in Qg excepted those belonging to 
the interface between Qg and its counterpart. This interface is denoted by J R. As 
explained in [94], the test functions are null on the interface J g, as if Dirichlet bound- 
ary conditions were imposed to the RID. On this interface, displacements follow the 
shape of the modes w, according to Eq. (2.19). The hyper-reduction method gives 
access to reduced coordinates (7%),—1,..., that fulfill the following balance equations, 
for contactless problems: 


50339) 


GH) = $ yE) ve, YX E€ Qe (2.23) 
k=1 
f e(Wrj) : oe) dQ (2.24) 
Opr 
-f vej fada- f Prj Tu. NdS=0. (2.25) 
OR IRRNAINy 
Vj =l,...,n (2.26) 


The matrix form of the hyper-reduced balance equations reads: find y € R” such 
that 


20 2 Learning Projection-Based Reduced-Order Models 


N 
Gx) = J 9, (x) Gi, YX E Qe, (2.27) 
i=l 
@= Vy, (2.28) 
FUR (y) := VIP, 17 FV IP], (2.29) 
FUR (y) = 0, (2.30) 


where V[P, :] denotes a row restriction of matrix V to indices in P. The reduced 
Newton-Raphson step reads: 


N 
a(x) = VaN, Wee QR, (2.31) 
i=1 
F = V yt, (2.32) 
K#? := vip, :]" a TDIP, :] V, (2.33) 
K”E (yt — y) = -VIP, 17 FADIA], (2.34) 


where the reduced stiffness matrix K”* is computed by using solely the elements of 
the RID Qg. We assume that the matrix K”? is full rank. This assumption is always 
checked in numerical solutions of hyper-reduced equations. Rank deficiency may 
appear when the RID construction do not account for the contribution of a reduced 
basis dedicated to stresses. 

Once the RID is represented as a finite element mesh, this hyper-reduction scheme 
is intrusive solely for the linear solver involved in the Newton-Raphson step and its 
related convergence criterion. Nevertheless, the mesh of the RID has to include labels 
for the set P or its counterpart Z g. This counterpart is the set of degrees of freedom 
connected to elements of the original mesh that are not in the reduced mesh. 


Remarks: 


Here the most complex operations are indeed the computation of K”* and the 
solution of the reduced linear system of equations. They respectively scale linearly 
with card(P) n? and n>. Hence n? has to be small enough compared to N if we 
consider the computational complexity for the solution of sparse linear systems in 
the finite element method. 

Because of the spreading nature of interpolation points, most of the time, the RID 
is not a compact subdomain. 

The hyper-reduced order model is a kind of submodel where the displacements at 
the interface J g follow the shape of the modes w, according to Eq. (2.19). 

Finite element corrections for displacements and stresses can be easily computed 
over the RID once the reduced prediction have been achieved. This scheme is 
termed Hybrid Hyper-Reduction in [46]. 

A parallel programming of the hyper-reduction method has been proposed in [95]. 
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e Reduced order models not only save computational time, they save computational 
resources including energy consumption savings as explained in [90] and memory 
footprint [46]. 


Property 1: In linear elasticity, if K”* is full rank, the hyper-reduced balance 
equations are equivalent to an oblique projection of the finite element prediction 
qe R”: 


I’ := V[P, :]" K[P, :], (2.35) 
q = V T VNI q, (2.36) 
and II’ = NF q, (2.37) 


with K q = F. Hence the hyper-reduced prediction of the reduced vector y is a 
minimizer for f (f): 


y* eR", f(y) =|" (Vy* —q) l$. (2.38) 


Here H is a projector for elastic stresses in Qg according to the reduced test functions: 


N 
X Nik Vy -@); =f efri) : (6(G) —o(q)) dQ. (2.39) 


i=l 


The proof is straightforward. Here, K¥? = I" V. The Jacobian matrix for f 
reads J = V” IM’ V = (K”*)" K”®*_ If K”? is full rank, then J is symmetric 
definite positive and J -1 = (K#£)-! (K#£)-T. Then, both the minimization prob- 
lem and the hyper-reduced equation have a unique solution. The solution of the 
minimization problem is: 


gf = VJ V" ON" q, (2.40) 
= V (K¥8)! N7 q, (2.41) 
= V (K#8)! VIP, :]" K[P, :]q, (2.42) 
=V(K"*)"! VIP, :]" FIFI, (2.43) 
= @. (2.44) 


As an intermediate result, Eq. (2.42) is the oblique projection. 
In linear elasticity, the Céa’s lemma holds. Let us denote q° € R” the minimizer 
of the upper bound in Eq. (2.3) related to this lemma: 
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N 
° = in || — ai |i. 2.45 
q° = arg min, | 29:4 I (2.45) 
N 
P = > 9:9). (2.46) 


The best projection of the minimizer q° in the approximation space is denoted by 
Y p: 


y?” = argmin,-g.(q° — V g)'M(q° — Vg), (2.47) 
N 

P=) pV yP. (2.48) 
i=1 


Let us introduce an ideal reduced basis V° e R‘*” (It assumes that n is an 
ideal reduced dimension) such that: q° = V° y°, and V°TMV? = I, where M; j= 
(9;,9;)- Hence y” = VMV” y°. 

Property 2: In linear elasticity, the upper bound of approximation error is 
increased by a Chordal distance [101] between V and the ideal reduced basis V°: 


jaw — T° || < IX- TP I + lly l2 d” (V°, V), (2.49) 


where d°” (V°, V) is the Chordal distance between V° and V. 

Hence, the smaller the Chordal distance between the sub-space spanned by V and 
V°, the better the reduced prediction by using a Galerkin projection (When the RID 
covers the full domain). A certification of the reduced projection can be achieved, 
when all errors admit an upper bound, by following the constitutive relation error 
proposed in [47, 55]. 

The Chordal distance uses the principal angles 0 € R”, & € [0, 7/2[ for k = 


1,...,, computed via a full singular value decomposition: 
V” MV’ =U cos(@)U", UU =UV =], (2.50) 
d™ (V°, V) = |isin(®)||F, (2.51) 
Ue =n, (2.52) 
where || - ||7z is the Frobenius norm. Here, cos(@) and sin(@) are cosine and sine 


diagonal matrices. In addition the following property holds when a full SVD is 
computed: 
UU" =I, UU’ =I. (2.53) 


The proof of the previous property is straightforward by using the triangular 
inequality. We just need to prove that: 


0° — BP I] < ya d” (V°, V). 
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Hence, the proof is the following: 


IDe _ DP |? = yT (V° _ Vv Mv’)? M (V° = VV' MV’) y° 


= y” (L — VMV VMV’) y”. (2.54) 
Therefore 
ID — XP ||? = y? (A — UP cos(0)? U7) y° (2.55) 
= yp? U°(I — cos(0)°) Uy? (2.56) 
= y” U°sin(0)? U°T y° (2.57) 
= |isin(0) U” y° |2. (2.58) 


For all matrices A € R”*” and B € R”*” the following property holds: 
IABIiF < ||Alle IIBIF, 


and for a € R”: |lal| = |lall2. 
Thus: 
© — X? < lsin(0) h4 WU? yeli < sinh yeni. (2.59) 


Property 3: When the identity matrix is substituted for K in Eqs. (2.35) and (2.36) 
is known as the Gappy POD reconstruction [35] of truncated variables q[P]. The 
reconstructed vector in RY is: 


q = V (VIP, 31’ VIP, :D VIP, 31” qiP]. (2.60) 


This Gappy POD reconstruction is useless for displacement variables because the 
oblique projection in Eq. (2.36) is a direct outcome of the hyper-reduced prediction. 
But such a reconstruction is very convenient for stress variables that the hyper- 
reduced scheme forecasts only on Qg. The reconstructed stress variables reads: 


X = V° (WP? , 1 VIP, Do WTP’, I? gp], (2.61) 


where P” is the set of all stress indices available in Q r- Since the RID contains 
interpolation points for V”, these points are included in P” , therefore the truncated 
matrix V” P, :] is full column rank and the reconstruction is a well posed problem. 
Remark about the RID construction and the DEIM: If the RID contains solely 
the elements connected to interpolation points related to the reduced basis V, such 
that P = P”, then the Gappy POD gives the interpolation scheme of the DEIM: 


FPEM = V (VIP", D! qP“). (2.62) 


But, when considering the hyper-reduction scheme, one can observe overfitting in 
the sense that the train set of displacement is very well approximated by the DEIM 
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reconstruction, but the hyper-reduced predictions are not accurate. For this reason, we 
recommend the use of the additional reduced basis V7” and the related interpolation 
points. 

Various applications of the hyper-reduction method using a RID have been devel- 
oped for: 


thermal problems in structures or solids, in [88], 

boundary element models [93], 

reduced simulations of sintering processes [99], 

ductile damage predictions, including unstable localisation of strains [97], 
reduction of multidimensional domains, when space variables are an Euclidean 
space of arbitrary dimension D > 3 [98], 

simulation of viscoelastic-viscoplastic composites materials [74], 

model calibration in plasticity of materials [37, 46, 96], 

contact problems using Lagrange multipliers [38, 62], 

arc length algorithm for buckling problems or strain localisation [59], 
micromorphic continua including higher order stress fields [48]. 


2.3.5 Hyper-Reduction via Empirical Cubature 


To assemble the linearized equations of the reduced Newton algorithm (2.13) when 
using the ROM in the online phase, hyper-reduction techniques via empirical cuba- 
ture aim to compute the costly integrals over the high-fidelity domain by replacing 
the high-dimensional quadrature formula by a low-dimensional reduced quadrature 
with positive weights. The ECSW [36], ECM [45] and LPEQP [100] are methods 
implementing such reduced quadratures. In this section, we present the ECM, more 
details are available in [19]. 

We consider the high-fidelity model described in Sect. 2.2. The integrals involved 
in the assembling of the linearized Eq. (2.7) make use of high-fidelity quadrature 
formulas. Apply such quadrature to the reduced internal forces vector: 


PGs fo (€(@), y) (x, t) : € (Wi) a) 
ne (2.63) 
=J} oo (eÔ), y) ar t) € (Wi) a,l <i <n, 


ecE k=1 


where E denotes the set of elements of the mesh, ne the number of quadrature points 
for the element e, œp and x; are the quadrature weights and points associated to e. 
The total number of quadrature points is denoted Ng. 

The ECM aims to approximate the high-fidelity quadrature by a reduced quadra- 
ture with positive weights, which, when applied to the reduced internal forces vector, 
writes 
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F(t) ax dye (eÔ), y) Gr, 1): € hi) Gv), 1 <i sn, (2.64) 
k'=1 


where ôy > 0 and X, are respectively the reduced quadrature weights and points, 
and ng < NG is the length of the reduced quadrature. 

Denote fy := o (€(u(q/jn)41)s Y) : € (Warni) 1 <q < nN.. where // and % 
are the quotient and the remainder of the Euclidean division. Denote as well Z”¢ 
a subset of [1; Ng] of size ng and Jz» € R"^e*ne and b e N’% such that for all 
1 <q <nN,.andall 1 <q' < ng, 


Jz = Czo) , b= ( f fa) , (2.65) 
4 l<q<nN-, q’€Z"G Q l<q<nN. 


where a denotes the q'-th element of Z”°. We remind that n is the number of 
NG 
POD modes, see Sect. 2.3.3.2. Let ô € R*G, (Jzro ô), = 5 qo (€(Uq@sjny+1) y) 
q'=l 
(xz"c ) z€ (Waan) Ezy) l<qŚ<nN., is a candidate approximation for 
q 


f o (elun) y) zE (Waan) = b4 , 1 < q < nN.. The problem of finding 
2 


the more accurate reduced quadrature formula of length ng for the reduced internal 
forces vector is: 


: (2.66) 


@,Z"°) = ar min Jznc@! — b 
( ) = arg y gr E y HZ" 
where ||- || denotes the Euclidean norm. Minimizing the length of the reduced quadra- 


ture formula as well leads to a NP-hard problem, which solution can be approximated 
using a Nonnegative Orthogonal Matching Pursuit algorithm, see Algorithm 2.2. 


Algorithm 2.2: Nonnegative Orthogonal Matching Pursuit. 


Input : J, b, tolerance €Op.comp.; Xk, | < k < NG 
Output: Ôk, Xks 1l<k<d. 


1 Set Z = Ø, k' = 0, ô = O and ro =b ; // initialization 
2 while ||rxll2 > €Op.comp. llbll2 do 
3 Z < Z U max index CERN rw) 
4 | ô< arg min |b — Jzô' |l 
ô'>0 
5 ry+ı < b — Jzô 
6 k <k +1 
7 end 
8d <k 
9 XK =xz,,l<k<d 
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In Algorithm 2.2, Jii; nc] satisfies the definition (2.65) with Z"° = [1; Ng]. The 
positivity of the weights of the reduced quadrature preserves the spectral properties 
of the operator associated with the high-fidelity problem, see [19, Remark 1]. 


2.3.6 Computational Complexity 


In this section we restrict our attention to elliptic problems or to linearized problems. 
The bilinear part of the weak form for finite-dimensional solution spaces is a matrix. 
When using Finite Element solution space, this matrix is sparse. But when using a 
reduced solution space, this matrix is usually a full matrix. Therefore, computational 
complexity of the finite element prediction is the complexity of the solution of sparse 
linear system. It scales linearly with N w*, where w is the band width of the sparse 
matrix. For the reduced prediction, the solution of a full linear system scales linearly 
with n>. We recommend to restrict linear model reduction schemes, with or without 
hyper-reduction, to reduced dimension n lower than N!/°, otherwise the solution of 
reduced equation will have a computational complexity similar to the complexity of 
the finite element model. This recommendation does not concern explicit solvers. 


2.4 Nonlinear Manifold Learning for Projection-Based 
Reduced-Order Modeling 


Consider a parametrized variability, and a set of snapshots generated using the high- 
fidelity model over a sampling of the parameter domain. The parametrized problem is 
said nonreducible when applying a linear data compression over this set of snapshots 
leads to a ROB containing too many vectors for the online problem to feature an 
interesting speedup. Formally, this happens when the Kolmogorov n-width d,,(M) 
decreases too slowly with respect to n, where we recall that n is the cardinality of 
the ROB, 


d,(M):= inf sup inf ||u — vll, 2.67 
(M) ei ee vile (2.67) 


with the Grassmannian Gr(n, H) being the set of all n-dimensional subspaces of 
H and H, € Gr(n, H) the subspace spanned by the considered ROB. Qualitatively, 
the solution manifold M covers too many independent directions to be embedded 
in a low-dimensional subspace. To address this issue, several techniques have been 
developed: 


e Problem-specific methods tackle the difficulties of some specific physics problems 
that are known to be nonreducible, such as advection-dominated problems which 
have been largely investigated, for instance in [16, 49, 85]. 
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Online-adaptive model reduction methods update the ROM in the exploitation 
phase by collecting new information online as explained in [102], in order to limit 
extrapolation errors when solving the parametrized governing equations in aregion 
of the parameter space that was not explored in the training phase. The ROM can 
be updated for example by querying the high-fidelity model when necessary for 
basis enrichment [18, 44, 56, 80, 88]. 

ROM interpolation methods [6—9, 24, 64-68, 75, 76] use interpolation techniques 
on Grassmann manifolds or matrix manifolds to adapt the ROM to the parameters 
considered in the exploitation phase by interpolating between two precomputed 
ROMs. 

Dictionaries of basis vector candidates enable building a parameter-adapted ROM 
in the exploitation phase by selecting a few basis vectors. This technique is pre- 
sented in [54, 72] for the Reduced Basis method. 

Nonlinear manifold ROM methods [57, 63] learn a nonlinear embedding and 
project the governing equations onto the corresponding approximation manifold, 
by means of a nonlinear function mapping a low-dimensional latent space to the 
solution space. 

Dictionaries of ROMs rely on the construction of several local ROMs adapted to 
different regions of the solution manifold. These local ROMs can be obtained by 
partitioning the time interval [32, 33], the parameter space [33, 34, 42, 44, 51, 
52, 82], or the solution space [10, 11, 27, 40, 77, 82, 92]. 


In the following Sects. 2.4.1 and 2.4.2, we provide more details on the last two 
entries of the previous list. 


2.4.1 Nonlinear Dimensionality Reduction via Auto-Encoder 


Nonlinear manifold learning means that the solution manifold is approximated by 
a domain in the ambient solution space that is not included in a low-dimensional 
vector subspace, as illustrated in Fig. 2.3. 

Let us consider a formal representation of parabolic and nonlinear Partial Dif- 
ferential Equations (PDEs) that are parameterized with respect to some physical 
parameters of interest. We mean by physical parameters the parameters that appear 
directly within the equations such as the boundary conditions, the viscosity for fluid 
mechanics (henceforth the Reynolds number), the time step for dynamical systems 
of fluid flows or infectious diseases, etc. These parameters are denoted u without any 
loss of generality as introduced in the preceding section. The formal representation 
of the equations is given as follows: 


Ou 
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Nonlinear data compression 


Variability with large influence on u Large solution manifold 


Parameter space Solution space H c L?(Q) 


Fig. 2.3 Nonlinear manifold learning 


Reduced order modeling based on nonlinear data compression techniques might 
be a solution for example when the described physical fields by the model equations 
require a large number of vectors in the ROM as specified above. Nevertheless there 
are cases where even if the physical solution fields are completely reducible, the 
Galerkin projection may not be appropriate for the model equations describing this 
physics. 

Convection-diffusion PDEs have this stability issue even more generally when 
considering Galerkin projection using finite element basis functions. In the literature, 
it is proved that the coherent structures of a turbulent, unsteady and in-compressible 
fluid flow are reproducible by a small number of POD basis functions. However, 
if these functions are used to solve the Newton-Raphson problem of the associated 
reduced order Galerkin dynamical system, then an instability appears as a function of 
the time. In the literature, many solutions are proposed to tackle this difficulty while 
keeping the reduced order approximation in a linear space spanned by the POD basis 
functions. We can refer to the Petrov-Galerkin technique, the least square minimiza- 
tion of the equations residual, the variational finite element method, etc. Recently, 
nonlinear approximations of the solution fields in a manifold of reduced dimension 
start to gain importance in the literature. In this case, the reduced order model is said 
to be a nonlinear projection based reduced model. Some authors introduced nonlin- 
ear approximations using Deep Learning approaches for projection based reduced 
models. We find in the literature more classical nonlinear approximations based on 
the Kernel POD technique. 

In what follows, we make a focus from the literature on Deep Learning projection 
based reduced models and their different possible formulations. More precisely, we 
are talking about Deep AutoEncoders (DAEs) from the domain of Deep Learning. 
DAEs are artificial neural networks formed of layers of spatial convolutions, nonlin- 
ear activation functions and linear systems called fully connected functions. These 
architectures are used to perform nonlinear dimensional reduction, following unsu- 
pervised data compression. Henceforth, the DAE allows to determine latent features 
within a set of given inputs. 
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We denote by h and g respectively the encoder mapping and the decoder mapping 
of a DAE. In general g is the transpose mapping of h. We denote by @ the reduced 
latent features inferred by h. The dimension of the latent features is equal to the 
intrinsic dimensionality of the manifold as stated in Remark 2.1 in [63]. This intrinsic 
dimensionality is in the current case the dimension N, of the vector of parameters 
u, Which may include the time variable also. 

We note that the reduced latent features of a DAE are not parameterized variables 
in general. In other words they can be seen as non-parametric features associated with 
a given set of inputs. This formulation is interesting in the framework of projection- 
based model reduction, where the associated parameters are given straightforwardly 
by the physical equations. Hence, knowing the variable parameters within the inputs 
data helps only with the determination of the intrinsic dimension of the manifold. 

In the literature, we find two different formulations for DAEs projection based 
reduced models. 

The first formulation was proposed by Kashima [53] and Hartman et al. [43]. 
It is very analogical to the Galerkin formulation of projection based reduced mod- 
els: given h and g such that goh:u —> t and T ~ i, then the reduced model is 
formulated as follows: 


ð 
zem =ho fog@wu)) , a(t =0) = ht = 0)) (2.69) 
=h(f(g(@(H)), n). (2.70) 


In the above formulation, the authors relied on the following three points in order to 
set the time derivative of the latent features equal to the right hand side of Eq. (2.69). 


a 
® 978 CWD belongs to the manifold described by u —> ù(u), 


a ð = 
eh (Zecu) = Th (e@(u))), 
ehog= Ty, 


Remark 2.1 The first two above items are hypothesis that are fulfilled in the case 
where h and g are linear or affine functions. The last item is fulfilled theoretically by 
the inputs data compression using parameters optimisation of the DAE architecture. 


The second formulation of DAEs projection based reduced models is proposed 
in [63], where a least square minimisation of the residual of parabolic PDEs because 
of the decoder approximation is performed. Then, the reduced model is formulated 
as follows in order to determine the reduced latent features: 


ð 
FEH = argming cx |I CUID — FEEW), WI, (2.71) 


where (u) is the time derivative of @(u), ||.||, denotes the mean square norm or the 
euclidean norm and, J is the Jacobian matrix of the decoder mapping which belong 


30 2 Learning Projection-Based Reduced-Order Models 


element-wise to the tangent space to the solutions manifold at a given point. J is 
expressed as follows: 
J:@>Ve@. 


In this second formulation, the authors do not suppose that the velocity of the decoder 
approximation is in the manifold of the solutions because mathematically it belongs 
to the tangent space to the manifold at a given point. Hence, they claim that encoding 
the decoder approximation will produce a poor approximation by the reduced model. 


2.4.2 Piecewise Linear Dimensionality Reduction 
via Dictionary-Based ROM-Nets 


Parts of this section has been inspired from the authors previous work [30]. 
Piecewise linear manifold learning means that the solution manifold is approxi- 
mated by a dictionary of local linear subspaces, as illustrated in Fig.2.4, where we 
denote M the solution manifold. 
The solution manifold is partitioned to get a collection of subsets My C M that 
can be covered by a dictionary of low-dimensional subspaces, enabling the use of 
local linear ROMs. If {Mg}xequ; xq is a partition of M, then: 


vk € [l; K], VN €e N*, dy(My) < dy(M). (2.72) 


The concepts of ROM-net and dictionary-based ROM-net are introduced in [27], 
which we present in this section. Suppose we dispose of an already computed dic- 
tionary of ROMs for the parametrized problem (2.4), where each element of the 
dictionary is a ROM that can approximate the problem on a subset of the solution 
manifold M. A dictionary-based ROM-net is a machine learning algorithm trained 
to assign the parameter jz € X to the ROM of the dictionary leading to the most accu- 


Piecewise linear data compression 


Variability with large influence on u Large solution manifold 


Parameter space Solution space H c L? (Q) 


Fig. 2.4 Piecewise linear manifold learning 
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K d 


Fig. 2.5 Exploitation phase of a dictionary-based ROM-net. K local ROMs, combined with a 
classifier Cx for automatic ROM Cx (u) recommendation, are used to predict the quantity of 
interest Z (u) 


rate reduced prediction. This assignment, called model recommendation in [77], is 
a classification task, see Fig. 2.5. 

The dictionary of ROMs is constructed in a clustering stage, during which snap- 
shots are regrouped depending on their respective proximity on M, in the sense of 
a particular dissimilarity measure we introduced in [29] and [26]. The dissimilar- 
ity between two parameter values 4, u’ € X, denoted by ô(u, u’), involves the sine 
of the principal angles between subspaces associated to the solutions of the HFM 
u(u), u(r’) E€ M, see [26, Definition 4.10]. Applying a k-medoids clustering algo- 
rithm on the solution manifold M using the dissimilarity ô leads to an optimal parti- 
tioning for a dictionary of local ROMs, in a sense introduced in [26, Property 4.13]. 
We refer to the remaining of [26] for the description of a practical efficiency criterion 
of the dictionary-based ROM-net, which enables to decide, before the computation- 
ally costly steps of the workflow, if a dictionary of ROMs is preferable to one global 
ROM, and how to calibrate the various hyperparameters of the ROM-net. 


Remark 2.2 Importance of the classification. One could argue that the classification 
step can be replaced by choosing the cluster k for which the dissimilarity measure 
ô(u, Ak) between the parameter u and the cluster medoid jz; is the smallest. However, 
we recall that the computation of the dissimilarity measure requires solving the 
HFM at the parameter value u, which would render the complete model reduction 
framework useless. Hence, the classification step enables to bypass this HFM solve 
and directly recommend the appropriate local ROM. 

As briefly mentioned in the introduction of Sect.2.4, local ROMs can be con- 
structed by partitioning the parameter space [33, 34], in which case the classification 
step is not required: the cluster affectation is made by computing distances directly 
in the parameter space. In other cases, partitioning in the solution space can be 
considered without requiring a classification step [10]. Consider a time-dependent 
problems where the initial condition is not a parameter of the problem, and sup- 
pose an efficient computation of the clustering distance in the solution space based 
on the reduced solution at the previous time-step. Then, local basis affectation and 
switching is possible without requiring classification. 


The training of the classifier can be difficult when working with physical fields: 
simulations are costly, data are in high dimension and classical data augmentation 
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techniques for images cannot be applied. Hence, we can consider replacing the HFM 
by an intermediate-fidelity solver for generating the data needed for the training of 
the classifier, by considering coarser meshes and fewer time steps. We point out that 
the HFM should be used at the end for generating the data required in the training of 
the local ROMs. We propose in [28] improvements for the training of the classifier 
in our context by developing a fast variant of the mRMR [83] feature selection 
algorithm, and new class-conserving transformations of our data, acting like a data 
augmentation procedure. 


2.5 Iterative and Greedy Strategies 


For the sake of the presentation, we have separated the offline and online phases, 
where the Reduced-Order Model is learnt, then exploited. Actually, more involved 
strategies exist, where the ROM is constructed in a iterative fashion. The Reduced 
Basis Method [84] is a greedy method, where the ROB is constructed by a sin- 
gle snapshot, corresponding to a randomly chosen parameter value, and the ROB 
is enriched by the parameter value that maximizes the error made by the current 
ROM. In complex parameter dependencies, the hyper-reduction scheme can be sim- 
ulateously constructed as the ROB grows, see [31] for such a scheme, with the EIM 
as hyper-reduction. This greedy construction as been extended to time-dependant 
problems in the POD-greedy method [41], and simultaneous hyper-reduction con- 
struction have also been proposed [88]. 

Such iterative strategies rely on a efficient computation of the error made by the 
ROM. Error estimation is investigated in the next chapter. 
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Chapter 3 A) 
Error Estimation ga 


3.1 Confidence and Trust in Model-Based Engineering 
Assisted by AI 


Consider first data-based machine learning techniques. They rely on large sets of 
examples provided during the training stage and do not learn with equations. Dealing 
with a situation that do not belong to the training set variability, namely an out- 
of-distribution sample, can be very challenging for these techniques. Trusting them 
could imply being able to guarantee that the training set covers the operational domain 
of the system to be trained. Besides, data-based AI can lack in robustness: examples 
have been given of adversarial attacks in which a classifier was tricked to infer a 
wrong class only by changing a very small percentage of the pixels of the input 
image. These models often also lack explainability: it is hard to understand what is 
exactly learned, what phenomenon occurs through the layers of a neural network. 
In some cases, information on the background of a picture is used by the network 
in the prediction of the class of an object, or bias present in the training data will be 
learned by the AI model, like gender bias in recruitment processes. Addressing these 
limitations is the subject of the Program Confiance.ai,! regrouping French academic 
as well as industrial partners from defense, transports, manufacturing and energy 
sectors. 

Model-based engineering enjoys better explainability—since the reference-model 
equations are known and used, and robustness—when the reference-model is well- 
posed. Concerning trust, in our projection-based reduced-order modeling context, the 
prediction is in general deterministic, and strict error bounds can be derived in par- 
ticular cases. This is a main difference with AI-based models, which use notions like 
confidence intervals and predictive uncertainties, expressed as probability results. In 


 https://www.confiance.ai/. 
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the remainder of this chapter, we present error bounds and indicators in projection- 
based reduced-order modeling, depending on the complexity of the underlying equa- 
tions. 


3.2 In Linear Elasticity and for Linear Problems 


Parts of this section has been inspired from the authors previous works [7] and [9]. 
We suppose that the problem of interest has the following discrete variational form, 
depending on a parameter u in a parameter set P: for a finite-dimensional space V 
of dimension N (with N > | resulting, e.g., from discretization), find u, € V such 
that 
En : ap (tu, V) = b(v), Wwe y, (3.1) 


where a,, is an inf-sup stable bounded sesquilinear form on V x V and b is a 
continuous linear form on V. We define the Riesz isomorphism J from V’ to V 
such that for all / € V’ and all u € V, (JL, u)qy = 1(u), where (-, -)y denotes the 
inner product of V with associated norm || - |y and V’ the dual of V. We denote 
~ ja, (u, v)| 
Bu i= ce 
ueVyery |lullvllullv 
itive lower bound of £,,. For simplicity, we consider that the linear form b is inde- 
pendent of the parameter u. The extension to -dependent b is straightforward. 
Applying the Galerkin method on the linear problem (3.1), using a ROB 
(Wi)i<i<n E€ R’** as done in Sect. 2.3.2 leads to finding u,, € V, such that 


> 0 the inf-sup constant of a, and Bu a computable pos- 


A 


Ê, : ap ûn, uj) = blu;), Vie {l,....n}. (3.2) 


The approximate solution on the ROB is written as (2.16): 
fin = Dov Wi. (3.3) 
i=l 


We assume that the sesquilinear form a, depends on p in an affine way, namely 
there exist d functions a; (4) : P — C and d y-independent sesquilinear forms a, 
bounded on V x ‘V such that 


d 
a (u,v) = Yaga (u,v), Yu,ve Yy, (3.4) 
k=1 


which enables the ROM to be online-efficient. 
Under the current assumptions, the following error bound holds (see 
[16, Sect.4.3.2]): for all u E€ P, 
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luy — alley < E) = Be WGyaylly, 
Ñ d 
> (3.5) 
=B,'|Got >) J aty{ Gui]. 


i=1 k=1 y 


with G, the linear map from Y to VY such that V 5 u œ> G,u := J (ay (u, ‘)- b) € 
V; G, inheriting the affine dependence of a, on u since, for all u € VY, 


d d 
Guu = —Jb +) o% Jak(u, +) = Goo + of Gxu, (3.6) 
k=1 k=1 


where Goo := —Jb e V and Gyu := Jag(u, -) € V forall k € {1,..., d}. 
The error bound (3.5) can rewritten in an equivalent way as 


Ñ d 
ĉu) = Ba (Goo, Goo)v + 2Re » > y/o (Gui, Goody 
i=l k=l 


1 


Ñ ad ? (3.7) 
+ 5 >» v agy" a" (Gru, Guj)y |, 
i, j=1k,l=1 
z 1 
= 3’ (8 + 2Re(s' iu) +H Su), 
where 5 := ||Goollv, s and x, are vectors in C% with components s; := (Gti, 


Goo)y and (£,)7 := ağ y;", and S is a matrix in C’" with coefficients S; y := 
(Gzui, Giuj)y (with I and J re-indexing respectively (k, i) and (/, j), for all 1 < 
k,l < dandall1 < i, j < n). The t superscript denotes the transposition. The vector 
s and the matrix S depend on the ROB (W;)ı<i<n but are independent of u, hence can 
be precomputed; while the vector x,, depends on the reduced basis approximation 
ii, Via the coefficients y“. A lower bound Ên of the stability constant of a, is also 
computed in complexity independent of N (which is possible, for example, by the 
Successive Constraint Method, see [10, 13]). 

We would like to stress that &; (u) = &2() (in infinite precision arithmetic): the 
indices 1 and 2 are used to denote two different ways to compute the same quantity. 
In particular, 6; (jz) is not online efficient, while & (jz) is. 


3.3 In Nonlinear Mechanics of Materials 


The remaining of this section is inspired from the authors work [8]. 
We look here for an efficient error indicator in the context of general nonlinearities 
and nonparametrized variabilities in nonlinear structural mechanics. The generality 
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of the assumptions and the complexity of the model lead us to search for quantities 
correlated to the error made by the ROM with respect to the HFM, instead of rigorous 
error bounds as considered in the previous section. 

The problem of interest is the same as described in Sect. 2.2, and we suppose that 
we have constructed a ROM following the methods described in Sects. 2.3 and 2.3.5, 
namely used POD for data compression and ECM for operator compression. 

The quantity of interest is the accumulated plastic strain over the complete struc- 
ture. Since this is a dual quantity, the ROM do not provide directly a prediction 
over the structure, but only at the reduced quadrature points selected by ECM, see 
Sect. 2.3.5. The Gappy-POD can be used for to recover the dual quantity of interest 
over the rest of the structure, see [8, Algorithms 3 and 4] for a presentation on the 
present context, and [11] for in seminal paper on Gappy-POD. 

A quantification for the prediction relative error of the accumulated plastic strain 
is defined as 


it ifipa #0, 


IPullr2@) 
E? := lpa-ĵ (3.8) 
x Pu-Pullp2 7 : 
H — ro otherwise, 
max || Py ll,2(e) 
HEPoff. 


where p,, and p,, are respectively the high-fidelity and reduced predictions for the 
accumulated plasticity field at the variability u, and Por, is the set of variabilities 
encountered during the offline stage. We underline the fact that p,, is the reduced 
prediction over the complete structure, hence after applying the Gappy-POD recon- 
struction. 

Define the ROM-Gappy-POD residual as 


IPu-Byll -e » 
IAE if|| P,,\l2 Æ 0, 


SP =) |p 51) P (3.9) 
H Pu Pie otherwise, 
max Ip, [2 
MEP off, 


where p, is the reduced prediction (after applying the Gappy-POD) taken at the 
reduced quadrature points (Pp k = Py (Xe), 1 < k < mP), Pu is the vector of the 
accumulated plastic strain as computed by the constitutive law solver at the reduced 
quadrature points during the online stage, and ||- ||2 denotes the Euclidean norm. 
Notice that in the general case, p,, # p,,: this discrepancy is at the base of our 
proposed error indicator. 

Notice that the relative error E% involves fields and L?-norms whereas the ROM- 
Gappy-POD residual &!, involves vectors of dual quantities in the set of reduced 
integration points and Euclidean norms. In (3.9), || p,, — P,.||2 is the error between the 
online evaluation of the accumulated plastic strain by the behavior law solver: p,,, and 
the reconstructed prediction at the reduced integration points X,: p,,, 1 < k < m”. It 
is explained in [8, Sect. 4.1] that || p,, — p,,|l2 is also the residual of the least-square 
optimization involved in the online stage of the Gappy-POD. 

Let B € R””*”” such that Bpi = YP Âk), 1 <k <m?,1 <i <n’, K :={py, 

for all possible variabilities u} and d(K, W)72(q) := aip irit. lu — wllzz(@), with 
vE 
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W a finite-dimensional subspace of L?(Q). The following propositions and corollary 
are proven in [8, Sect.4.1]. 


Proposition 3.1 There exist two positive constants Cı and Cz independent of u (but 
dependent on n? ) such that 


2 m PS 
12@) < CillBzy — Bylld + Cill Py — Bull + C2d(K, Spanth? hzizne) 2o) 
(3.10) 


| Pu — Pp 


Proposition 3.2 There exist two positive constants K, and K independent of u 
such that 


m A ~ 2 A 
Bu — Bull? < Ki | Pu — Palio + Kolpa PN- (3.11) 


Corollary 3.3.1 Suppose that the reduced solution is exact up to the considered time 
step at the online variability u: py, = Py, in L? (Q). In particular, the behavior law 
solver has been evaluated with the exact strain tensor and state variables at the inte- 
gration points xx, leading to Pu (Âk) = Pu(Xk), 1 < k < m’. From Proposition 3.2, 
By — Bulls =0, and &} = 0. 


We observe that in practice, the evaluations of the ROM-Gappy-POD residual 
El (3.9) and the error E A (3.8) are very correlated in our numerical simulations. The 
idea is to exploit this correlation by training a Gaussian process regressor for the 
function &/, +> EŻ. At the end of the offline stage, we propose to compute reduced 
predictions at variability values {2;}\<;<y, encountered during the data generation 
step, and the corresponding couples (El, Eh). 1 <i < N.. A Gaussian process 
regressor is trained on these values and we define an approximation function 


8 > Gpr” (8f), (3.12) 


for the error E at variability jz as the mean plus 3 times the standard deviation of the 
predictive distribution at the query point S/: this is our proposed error indicator. If the 
dispersion around the learning data is small for certain values &/,, then adding 3 times 
the standard deviation will not change very much the prediction, whereas for values 
with large dispersion of the learning data, this correction aims to provide an error 
indicator larger than the error. We use the GaussianProcessRegressor python class 
from scikit-learn [17]. Notice that although some operations in computational com- 
plexity dependent on N are carried-out, we are still in the offline stage, and they are 
much faster than the resolutions of the large size systems of nonlinear equations (2.8). 
If the offline stage is correctly carried-out and since &4 is highly correlated with the 
error, only small values for &/; are expected to be computed. Hence, in order to train 
the Gaussian process regressor correctly for larger values of the error, the reduced 


Newton algorithm (2.13) is solved with a large tolerance ey = 0.1. We call these 
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HF 


solutions 


offline 
variability 


data comp. modes and reduced 
reduced solver 
operator comp. integration 


data generator 
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u 
Hi i Bo 
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Pr varT co po Pu, 0, 
E7 > Gpr Ep Epa and error Puis Cui 


Fig. 3.1 Workflow for the offline stage with error indicator calibration [8] 


operations “calibration of the error indication”, see Algorithm 3 for a description 
and Fig. 3.1 for a presentation of the workflow featuring this calibration step. 


Algorithm 3: Calibration of the error indicator. 
Input : Outputs of the data generation, data compression and operator 
compression steps of Section 2.3.3 
Output: Approximation function &! > Gpr” (84) of the error 
1 Set Z = Ø, k' = 0, ô = O and rọ =b ; // initialization 
2 fori =1..., N; do 
3 | Construct and solve the reduced problem (2.13) with ROM = 0.1 
4 | Compute the reconstructed plasticity p,,, using Gappy-POD and &/,, 
5 | Compute the error E}, using (3.8) 
6 | X < XU (8, El) 
7 Construct an approximation function &/, > Gpr” (84) of the error Ef using a 
Gaussian process regression and the data from X 


We recall that in model order reduction, the original hypothesis is the existence 
of a low-dimensional vector space where an acceptable approximation of the high- 
fidelity solution lies. The hypothesis is formalized under a rate of decrease for the 
Kolmogorov n-width with respect to the dimension of this vector space. The same 
hypothesis is made when using the Gappy-POD to reconstruct the dual quantities, 
which are expressed as a linear combination of constructed modes. For both the 
primal and dual quantities, the modes are computed by searching some low-rank 
structure of the high-fidelity data. The coefficients of the linear combination for 
reconstructing the primal quantities are given by the solution of the reduced Newton 
algorithm (2.13). After convergence, the residual is small, even in cases where the 
reduced order model exhibits large errors with respect to the high-fidelity reference: 
this residual gives no information on the distance between the reduced solution and 
the high-fidelty finite element space. 
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However, in the online phase of the ROM-Gappy-POD reconstruction (see [8, 
Algorithm 4]), the coefficients P, (the accumulated plastic strain computed by the 
constitutive law solver during the online stage) contain information from the high- 
fidelity behavior law solver. Moreover, an overdetermined least-square is solved, 
which can provide a nonzero residual that implicitly contains this information from 
the high-fidelity behavior law solver: namely the distance between the prediction 
from the behavior law and the vector space spanned by the Gappy-POD modes 
(restricted to the reduced integration points): this is the term || Bz„ — p ull2 in (3.10). 
Hence, the ability of the online variability to be expressed on the Gappy-POD modes 
is monitored through the behavior law solver on the reduced integration points. When 
the ROM is solved for an online variability not included in the offline variabilities, 
then the new physical solution cannot be correctly interpolated using the POD and 
Gappy-POD modes: hence, the ROM-Gappy-residual becomes large. 

From Proposition 3.2, if ||Bz, — p,,ll2 = ||P, — Pull is large, then the global 
error | Pur Pu | LO) and/or the error at the reduced integration points x, is large, 


which makes ||Bz,, — p ull2 a good candidate for a enrichement criterion for the 
ROM. A limitation of the error indicator can occur if the online variability acti- 
vates strong nonlinearities on areas containing no point from the reduced integration 
scheme, namely through the term Cd (K, Span{ y? }i<icnr jg) in (3.10). 

We recall that the error indicator (3.12) is a regression of the function 6, + Ef. 
In the online phase, we only need to evaluate &!; and do not require any estimation 
for the other terms and constants appearing in Propositions 3.1 and 3.2. 

Equipped with an efficient error indicator, we are now able to assess the quality 
of the approximation made by the reduced order model in the online phase. If the 
error indicator is too large, an enrichment step occurs: the high-fidelity model is 
used to compute a new high-fidelity snapshot, which is used to update the POD and 
Gappy-POD basis, as well as the reduced integration schemes. Notice that for the 
enrichment steps to be computed, the displacement field and all the state variables of 
the previous time step need to be reconstructed on the complete mesh Q to provide 
the high-fidelity solver with the correct material state. The workflow for the online 
stage with enrichment is presented in Fig. 3.2. 

We refer to [8] for more details on this subject, and detailed numerical applications 
in nonlinear structural mechanics for this error indicator and its ability to enrich a 
ROM in the online stage. 

Notice that another (noncertified) indicator in nonlinear solid mechanics with 
internal variables has been proposed in [1], aiming to approximate the dual norm 
of residuals in the same fashion as in the linear case described in Sect.3.2. For 
such nonlinear case, rigorous error bounds are not obtained: a gappy-POD-based 
approximate representation of the stress tensor is used, and the inf-sup constant 
evaluation has been replaced by a normalization of the residual using the norm of 
the Riesz elements for the external loading. 
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Fig. 3.2 Workflow for the online stage with enrichment [8] 


3.4 In Computational Fluid Dynamics 


In the following section, we present a priori error estimates due to the POD-Galerkin 
approximation applied to fluid dynamics equations in particular and parabolic nonlin- 
ear PDEs in general. It is a theoretical result on the convergence of the POD-Galerkin 
reduced order model towards the high-fidelity semi-discretized equations in the sense 
of the spatial variable. The solution of these semi-discretized equations is denoted 
a” over a time interval [0, T] such that h = + and M is the cardinal of a Hilbert 
basis that is capable of generating the high-fidelity semi-discretized solution at some 
specified M time instants. This orthonormal basis can be obtained in an a posteriori 
fashion thanks to the snapshots POD method. It is denoted by (W;);=1,...,m- 

The following convergence result includes furthermore a discussion on the stabil- 
ity of the Galerkin projection technique applied to parabolic PDEs. This result has 
been developed and published in the following papers: [2-5]. 

In the literature, we can find many works around the problem of defining conver- 
gent and a priori upper bounds for the POD-Galerkin reduced models for parabolic 
equations. It is a subject of great interest, if we can quantify efficiently the error of a 
reduced solution # obtained with an approximation technique of the corresponding 
high fidelity semi-discretized solution 7”. This problem can be seen as a theoreti- 
cal confidence interval around a training point of an approximation model (details 
on some convergence results of the literature should be added). Let us denote by 
Q the open and bounded domain of the spatial variable, such that Q € R“, where 
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d=1 or 2.Letus consider the parabolic PDEs for which the weak formulation in 
the space of the solution u” spanned by the POD basis of cardinal M, is described 
as follows: 


e bisa trilinear form defined over [H!(Q)]¢ x [H'(Q)]4 x [H!(Q)]4, 

e a dissipating term defined as a bi-linear and coercive form a over [H 1(Q)]f x 
[H (QË, 

e a linear form F, defined over [L?(Q)]¢, 

e $ is the coercivity constant of the bilinear form a, 

e C, and C; are respectively the norms of a and b in the space [L7(Q)]¢, 


eK= | u’ leaozman 


The following result is proved: 


Theorem 3.1 [fn << M is the dimension of the truncated POD-Galerkin reduced 
solution @ that represent the training point (solution) u", then we can derive the 
following a priori upper bound of the [L?(Q)]¢— error between T(t) and Ñ” (t), Y 
t €[0, T]: : 

IE -DO lio < AM). (3.13) 


Where fi(n) is the remainder of the sum of a convergent series; fi(n) converges to 0 
when N converges to M. More precisely, f\(n) is a function of the remainder of the 
sum of the POD eigenvalues (Ax)x=1,...m Obtained from the Snapshots POD applied 
to M temporal snapshots of u": it is expressed into two different fashions: Either 


C2 x 
fi(n) = (1 +2C,K + 2) T Š he (3.14) 


€ k=n+1 


if (€ — 26+ 6C,K) < 0, fora strictly positive real number €, or 


2 M 
fin) = (r + (200K + <) T exp(T[2C,K + K°(1 + cop) > Àk, 
ST 3.15) 


if (€ — 2B + 6C, K) > 0, for all strictly positive real numbers €. 


The mathematical proof of Theorem 3.1 is based on the properties of the forms 
a and b, the application of the Young inequality and the Gronwall lemma. For the 
details of the proof, refer to [4]. 


Remark 3.1 The result of Theorem 3.1 is applicable in particular for the 1 D Burg- 
ers equation and the 2D unsteady and incompressible Navier-Stokes equations, by 
remarking the following: 


e When d = 1, H! (Q) C L® (Q): there exists (OR € R** such that Y v € H! (Q), 
lvlzæ < (OF IVvllzzo. 
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e When d = 2 ord = 3, H! (Q) C L*(Q): there exists cl e R*™ such that Vv € 
H'(Q), lvl < C4 WV vll 2). 
e When d = 2: there exists C € R** such that: 


1/2 1/2 
lullo < C lvllzzo Voll ce . 


e If we denote by S$ the square matrix of dimension M defined by: Sj; = 
(Wi, y jh HQ? then Y v in V” the solutions space spanned by the complete POD 
basis y we have the following inequality: ||| ~ig ye < VISI olli - Where 
|| S|| is the spectral norm of the matrix S. For details, refer to [14]. 


Remark 3.2 A result of stability with respect to time of the POD-Galerkin reduced 
model for nonlinear and dissipated parabolic PDEs can be derived from result 3.1. 
If denotes only the viscosity constant in this particular case (without any loss 
of generality), then the POD-Galerkin reduced model is stable with respect to time 
when u satisfies the following inequality: 


6C,K +e 
pe 
UZ 2B , 
where € is a strictly positive real number. 


Remark 3.3 More particularly, in the case of linear and dissipated parabolic PDEs, 
the stability condition of the POD-Galerkin reduced model becomes equivalent to 
€ 


u . The error of the POD-Galerkin reduced model with respect to the high 


2Cp 
fidelity training point can be estimated exactly as in the Céa lemma applied for 
elliptic and sesquilinear PDEs: if u > 1 then, 


Jat -DO < ieee T 5 à 
L(Y = 26 ms 


n=N+1 


Based on the same methodology, we propose in what follows an a priori error 
estimate and a convergence result for POD-Galerkin reduced model parameter wise. 
In other words, we show that when the parameters change with respect to the training 
ones, a confidence interval is obtained around the new test solution. The width of 
this confidence interval converges to the truncation error of the ROM at the training 
parameters. This parametric convergence result is formulated as follows: 


Theorem 3.2 Let us denote by i, a training solution associated with a train- 
ing parameter uo. So we suppose we have only one training point for the POD- 
Galerkin reduced model, without any loss of generality. If Y™ denotes the Hilbert 
basis obtained from the POD applied to the high fidelity training solution and ty, w 
denotes the truncated POD-Galerkin reduced solution that approaches the test point 
(solution) u’, in the reduced POD space of dimension n spanned by yy", then: 
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|t — p) OW < FEC) (1+ lu- uol”) + °C) lu- polt, (8-16) 


where f3” (n) is the remainder of the sum of a convergent series; f3” (n) converges 
to 0 when n converges to M. More precisely f3” (n) is a function of the remainder of 
the sum of the orthogonal projection coefficients of the characteristic function 1g, 
such that Q, tends to Q when x tends to 02: it is expressed as follows: 


M 


—1/2 
= B 
2 = (lo, -0l *+ Fewczan-n)] Ca [Vul DS (la. wey, 


k=n+1 
(3.17) 
where A and B are strictly positive constants of which the detailed expressions are 
given in [2]. 


The proof of the above theorem is published in [2]. It is based on the properties 
of the two forms a and b, the application of the Young inequality and the resolution 
of a nonlinear ordinary differential inequality of Ricatti type. 


3.5 A Note on Accuracy of a Posteriori Error Bounds 
and Round-Off Errors 


In this section, we explain why the online-efficient error bound (3.7) may be sensitive 
to round-off errors. 

In computers, real numbers are represented by a finite number of bits, called 
floating-point representation. Current architectures are optimized for the IEEE 754 
double-precision binary floating-point format. Let x and y be real numbers. When 
computing the operation x + y, the result returned by the computer can be different 
from its theoretical value. Whenever the difference is substantial, a loss of signifi- 
cance occurs. A well-known case of loss of significance is when x and y are almost 
opposite numbers. Suppose that x = —y. We denote by maxfl(x + y) the result that 
the computer returns when the maximal accumulation of round-off errors occurs 
when computing the summation. There holds 


|maxfl(x + y)| ~ 2e|x|, (3.18) 


where € is called the machine precision. In double precision, € = 5 x 107!6 (see [12, 
Sect. 1.2]). 

When implementing an algorithm, one should ensure that each step is free of such 
a loss of significance. In some cases, simply changing the order of the operations 
can prevent these situations. As an illustration, consider x = 1, y= 1+ 1077, and 
the operation x? — 2xy + y*. This is a sum of terms where the first intermediate 
result in the sum is 14 orders larger than the result. Therefore, a loss of significance 
is expected. The relative error of this computation is about 8 x 1074. Computing 
(x — y)*, which is the factorization of the considered operation, leads to a relative 
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error of about 10~°. Thus, the terms of the sum are only 7 orders larger than the 
results, leading to a less catastrophic loss of significance. In this specific case, the 
remedy consists in carrying out the sum before the multiplication. In our projection- 
based ROM context, the evaluation of the formula & suffers from such a loss of 
significance, as we now explain. 

We investigate the influence of round-off errors when computing the error bounds 
(3.5) and (3.7) for respectively 6; (u) and &2(jz). As observed in the previous para- 
graph, the computation of a polynomial using a factorized form is more accurate than 


p 2 
using the developed form, in particular at points close to its roots. Here, (B£: w) 


is a multivariate polynomial of degree 2 in X,, computed in a developed form, whereas 
the scalar product (G „uu, G pup )y used in the computation of &ı (u) is not devel- 
oped. The following holds (see [9, Proposition 2.2.1] for the proof) 


Proposition 3.3 Letu € P and let maxfl(6,,E, (u)), k = 1, 2, denote the evaluation 
of Bu Sx (u) when the maximum accumulation of round-off errors occurs. There holds 


maxfl(6,,81()) = 28e, 


x (3.19) 
maxfl(B,E2(u)) > 28€, 


where 5 = ||Gool|y and € is the machine precision. 


From this proposition, we notice that the online-efficient formula &2 (u) suffers 
from an important loss of significance. 

We present below an error estimator proposed in [9] that enjoys both accuracy and 
online-efficiency. Leto := 1 + 2dn + (dn). Fora given u € Pria and the resulting 
u, € Span{w, ..., Yn} solving the reduced problem (3.2), we define (w) e C” as 


the vector with components (1, £,,,,.%* , 2% %,,), where X,,, = a y” (we recall that 
H 


My? Hi 
y; are the coefficients of the reduced solution in the reduced basis, see (3.3), and 
ak the coefficients of the affine decomposition of a, in (3.4)), with 1 < J, J < dn 
(with J = i + n(k — 1) such that 1 < i < n,1 < k < d, and with J = j +n(l — 1) 
such that 1 < j < n,1 < I < d). We can write the right-hand side of (3.7) as a linear 
form in X (u) as follows: 


5°? + 2Re(s‘k,) + "Sk, = 2 i500); (3.20) 


p=1 


where t, is independent of u (as ô, s, and S are independent of jz) and x p(/4) is the 
p-th component of x (u). 

Consider the function of two variables (p, u) b> £, (u), for all p € {1,..., 0} 
and all u € P. We look for an approximation of this function in the form 


Yu EP, Yp € {1,0}, Xp(u) ~ J A WR plur), (3.21) 


r=1 
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for a certain parameter ô < o. The Empirical Interpolation Method (EIM) provides a 
numerical procedure to construct this approximation and to choose the interpolation 
points (see [6, 15]), which leads to the following formula for computing the error 
bound 


1 


2 


Eu) = Pa'h J GOV, | (3.22) 


r=1 


where V, = | Gy, tp, a and where Aĉ (u) and u, are provided by EIM, see [9, 
Sect. 3.2] for all the details of this derivation. There holds (see [9, Proposition 3.2.1]): 


Proposition 3.4 The computation of the formula E; is well defined, and this formula 
is online-efficient. 


Besides, &3 involving a linear combination of accurately computed scalar products 
(see (3.5)), it is not subject to the loss of significance encountered in &). 

We refer to [9] for more details on the notion of validity of a formula to compute an 
error bound, additional variants for accurate and efficient error bounds (including one 
featuring a stabilized EIM), as well as numerical illustrations for a one-dimensional 
linear diffusion problem and and three-dimensional acoustic scattering problem. The 
error bound &3 can be of particular interest in the following situations: (i) when the 
stability constant of the original problem is very small (this is the case in many prac- 
tical problems), (ii) when very accurate solutions are needed, (iii) when considering 
a nonlinear problem (for which, in some cases, no error bound is possible until a 
very tight tolerance is reached, see [18]). 
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Chapter 4 ®) 
Resources: Software and Tutorials EE 


4.1 Mordicus: Reduced-Order Methods Designed 
for Industrial Usage 


MORDICUS is the French acronym for “Méthodes d’Ordre RéDult Conçues pour 
des Usages induStriels”, translated to “reduced-order methods designed for industrial 
usage”. It is the name of a collaborative project that took place from 2018 to 2023, 
with the objective of developing a standard for a datamodel and basic computational 
treatment for reduced-order modeling in the French community. 


4.1.1 Mordicus Project and Consortium 


“MOR_DICUS” is a project of type FUI (Fonds Unique Interministériel), in which 
took part a consortium of ten partners of small and large companies, as well 
as research institutes and universities. More precisely, the partners are Électric- 
ité de France (EDF), Safran Group, ESI Group, Hexagon, Phimeca, Transvalor, 
CT Ingénierie, CEMOSIS, Armines and Laboratoire Jacques-Louis Lions (LJLL- 
UPMC), see Fig. 4.1. 

The founding is provided by BPI France, pôle systematic and pôle Alsace 
énergie Vie, see Fig. 4.2. 


4.1.2 Mordicus Library 


Mordicus consists in a Python library and a C++ application, where a datamodel and 
common basic algorithms and tools are proposed to numerically carry out surrogate 
modeling and projection-based ROM. 
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Fig. 4.1 Partners of the MOR_DICUS consortium 
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Fig. 4.2 Funders of the MOR_DICUS consortium 


The sources of the library are available on GitLab.com! and the user’s license 
is the GNU LGPLv3. A website containing an installation procedure, a description 
of the library, some tutorial and a documentation is also available.* The rest of the 
section is inspired from the “Numerical methods”? section of this documentation. 

The Mordicus library is constructed with the datamodel at its center, and the 
methods and algorithms act as functions that can modify the state of the datamodel. 
As detailed in Sect.2.3, Reduced-Order Model workflows involve an offline (or 
learning) stage, and an online (or exploitation) stage. The different steps can be 
vastly adapted and mixed together to produce any complex workflows, as long as 
the algorithm makes sense. For instance, a regressor can be trained directly on high- 
dimensional data, and the data compression step would not exist. Another example 
is the Reduced-Basis method, where the notions of offline and online stages merge 
together, since HFM and ROM resolutions are alternatively computed when con- 
structing the reduced-order basis. 

Mordicus provides a datamodel, standard interfaces and basic algorithms, as illus- 
trated in Fig.4.3. Based on Mordicus, applicative modules can then be developed 
and propose new workflows and algorithms. The library genericROM, described in 
Sect. 4.2, is an applicative module of Mordicus. 

Simple algorithms are proposed, and can be used by any applicative module. 


e SVD: for computing truncated singular value decomposition of lower triangular 
matrices. 


l https://gitlab.com/mor_dicus/mordicus. 
? https://mordicus.readthedocs.io. 
3 https://mordicus.readthedocs.io/en/latest/_methods/index.html. 
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Fig. 4.3 Organization and overview of the datamodel of the Mordicus library 


e ScikitLearnRegressor: contains examples for computing grid search cross 
validation and a customizable Gaussian Process Regressor using scikit-learn tools. 


e Interpolation: contains efficient time Interpolation tool used in the Mordicus 
datamodel. 


More details on the datamodel are provided in the next Section. 


4.1.3 Mordicus Datamodel 


The main feature of Mordicus is a datamodel adapted to reduced-order modeling. 
It has been constructed to facilitate collaboration, by proposing three main classes 
CollectionProblemData, ProblemData and Solution, supposed to be 
populated and handle in the same fashion by all users (see Fig. 4.3): 


e CollectionProblemData: The meta structure containing the complete data- 
model for a reduced-order model. 

e ProblemData: Containing a model for a physics problem: initial condition, 
loading, constitutive laws, solutions. 

e Solution: Containing the size, snapshots and reduced coordinated of solutions. 


These classes also contain numerous functions to iterate, modify and handle the 
data structure. We present a few important ones below: 


e CollectionProblemData functions: 


— DefineVariabilityAxes (): sets the axes of variability, can be strings 
for nonparametrized variability, or floats, 

— SnapshotsIterator (): returns an iterator over snapshots of solutions of 
a given name, 


56 4 Resources: Software and Tutorials 


— GetSnapshotsAtTimes (): returns an array containing all the snapshots of 
a given name, interpolated at a given time, 

— CompressSolutions(): compress the snapshots of solutions of a given 
name against the corresponding reduced-order basis, and update to correspond- 
ing solution.reducedCoordinates. 


Notice that the functions acting on solution objects automatically iterate over all 
the problemDatas included in the collectionProblemData. 
e ProblemData functions: 


— UncompressSolution(): uncompress the reducedCoordinates of a solu- 
tion of a given name, and update to corresponding solution.snapshots, 

— GetLoadingsOfType (): returns all loadings of a specific type in a list, 

— GetParameterAtTime(): returns the parameter value at a specitiy time 
(with time interpolation if needed). 


e Solution functions: 


— UncompressSnapshotAtTime (): uncompress the reducedCoordinates of 
the solution at a given time, and update to corresponding snapshots, 

— GetTimeSequenceFromSnapshots (): returns the time sequence from 
the snapshots dictionary, 

— ConvertReducedCoordinatesReducedOrderBasisAtTime (): 
converts the reducedSnapshot at a given time from the current reducedOrder- 
Basis to a newReducedOrderBasis using a projectedReducedOrderBasis. 


The data-model has also been thought to be agile and customizable, by allowing 
developers to propose other classes, in their applicative module, or variant classes of 
the ones contained in the subfolders of Containers (e.g. Loadings, Meshes). 


4.2 GenericROM Library 


The genericROM software consists in a Python library, acting as an applicative 
module of Mordicus, and developed at Safran. 

The sources of the library are available on GitLab.com* and the user’s license 
is BSD 3-Clause, a very permissive license that allows the users to use and redis- 
tribute the sources, with or without modification. A website containing an installation 
process, a description of the library, some tutorial and a documentation are also avail- 
able.’ 


4 https://gitlab.com/drti/genericrom. 
5 https://genericrom.readthedocs.io. 


4.2 GenericROM Library 57 


4.2.1 Main Available Methods 


We refer to Sect. 2.3 for an description of the main steps of projection-based reduced- 
order modeling methods. 

The simplest data compression method available in genericROM is the snapshot- 
POD, as described in Sect. 2.3.3.2. An implementation is available in parallel with 
distributed memory, by partitioning the domain Q in subdomains, as well as an 
incremental version. 

The only hyper-reduction method available to date in genericROM is the ECM, see 
Sect 2.3.5. A variant of th Nonnegative Orthogonal Matching Pursuit Algorithm 2.2 is 
implemented, where randomly chosen quadrature point can be added at each iteration 
to prevent the algorithm to fall into local minima. 

The physical problems available for reduction by genericROM are 


e nonlinear quasistatic structural mechanics, with possibibly very complex constitu- 
tive laws, centrifugal and temperature effects, as well as pressure and homogeneous 
Dirichlet boundary conditions; 

e transient thermal problems, with heat capacity and conductivity depending on the 
solution temperature in a nonlinear fashion, as well as convection heat flux and 
radiation boundary conditions. 


4.2.2 Noninstrusivity and Nonparametrized Variability 


In genericROM, projection-based ROMs can be constructed even when the snap- 
shots are generated by commercial software, as long as readers for the meshes and 
computed snapshots are available (or can be developed). The handling of projection- 
based ROM workflows without having to modify the assembling routines of the 
reference HFM code is the reason why the library is coined nonintrusive. To do so, 
the assembling routine of operators and right-hand sides are handling directly in 
genericROM, using the finite element engine of BasicTools® developped at Safran. 

This nonintrusive feature is rare in the literature, and has important advantages in 
an industrial context: expensive and already computed databases of snapshots can be 
used to construct a ROM. Besides, the ROM library can be updated without having to 
undergo complex certification of the evolution of reference HFM codes. Moreover, 
when these codes are commercial, intrusive implementations of projection-based 
ROM are not possible, or may require expensive contracts with code editors. An 
example of a ROM constructed with genericROM for a nonlinear transient thermal 
problem using snapshots computed by Abaqus is available in [2]. 

GenericROM can also deal with nonparametrized variability, in the sense that 
the variability at the origin of the differences in the snapshots computed in the data 


6 https://gitlab.com/drti/basic-tools. 
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generation step is not required to be formalized, or even known. This variability takes 
usually the form of parameters in the definition of the reference HFM, but can also 
be unknown when for instance the boundary conditions come from a first numerical 
simulation. Besides, one may want to compute a reduced model by imposing a 
variability outside a parametrization used to generate the database of snapshots. The 
nonparametrized variability feature comes with an efficiency trade-off, since, for a 
certain variability to be changed in a nonparametrized fashion, the ROM needs to 
assemble the uncompressed version of the corresponding terms in the weak form, 
before compressing them in the reduced problem (2.13). GenericROM offers the 
possibility to precompute all the parametrized variability, in order for the ROM 
to compute efficiently the corresponding terms in the reduced problem, with the 
consequence that such terms can only be updated following this parametrization. 
More details on the efficiency are given in the next section. The next two sections are 
inspired from the “Numerical methods” and “Tutorials”? section of the genericROM 
documentation. 


4.2.3 Precomputations for Efficiency 


In ROM workflows, the common measure of efficiency is the speedup, defined as 
the ratio between the computation duration of the HFM and the one of the ROM. For 
a given accuracy, the higher the speedup, the better. Various elements can be taken 
into account, including code optimization of the online stage. The methods enabling 
algorithmic complexity gains have been presented in Sect. 2.3.3. In this section, we 
give additional implementation elements and illustrations to give some intuition on 
the mechanisms at play. 

Depending on the considered equations and nature of variability, the goal is to 
precompute as many quantities as possible in the offline stage, to leave as few opera- 
tions as possible to the online stage. In the HFM, the assembling of the global tangent 
operator at each Newton iteration, using the high-fidelity quadrature, reads: 


DF a 
Pfu u -Xok ):K (et), Yu): €D] (xe); 1<i, j <N. (4.1) 


In the ROMs implemented in genericROM, the assembling of the reduced global 
tangent operator at each Newton iteration, using the reduced quadrature computed 
by ECM, reads: 


7 https://genericrom.readthedocs.io/en/latest/_methods/index.html. 
8 https://genericrom.readthedocs.io/en/latest/_tutorials/Mechanical/MecaSequential/index.html. 
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Fig. 4.4 HFM: a large 
sparse linear system is 


assembled and solved at each 
step of the Newton algorithm 
N = 
k 
ü se (uk) 
Fig. 4.5 ROM: a small n 
dense linear system is <> 
assembled and solved at each 
step of the reduced Newton n | || = | 
algorithm 
DF (pk ak+1 ak k 
Dù (âi) aa Fy (âi) 
Ng 
DF, F x A K ak : $ l<i.j< 
ai (a L N Xô [e (wi) : GUME Yu) ee (Wi) | (îe); <ij<n&N. 
g=! 
(4.2) 


The corresponding linear systems are illustrated on Figs. 4.4 and 4.5. In particular, 
n should be much smaller than N to obtain a significant speedup despite the fact that 
the linear system is sparse in the HFM and dense in the ROM, see Sect. 2.3.6. 

Based on Figs. 4.4 and 4.5, the speedup is clear for the resolution part of the linear 
system. We still need to illustrate how precomputations in the offline phase allow an 
efficient assembling of the reduced problem. Denote d the number of unknowns of 
second-order tensor dual quantities. In 3D, d = 6, for instance for the strain tensor, 
the unknowns are €11, €22, €33, €12, €23, €31. The reduced global tangent operator at 
each Newton iteration, using the reduced quadrature, can be written 


DF, G) x 

Di ij 

Ng d d 

we Ll [er (Wi) (%e)] X [Kim (EC), Yu) (%e)] [em (Wi) (%e)]. 1S ij <n. 
g=l l=1 m=1 


(4.3) 
In genericROM, the code for computing this quantity in the online stage is: 


reducedTangentMatrix = np.einsum(‘’g,lgj,glm,mgi->ij’ 
reducediIntegrationWeights, 
reducedEpsilonAtReducediIntegPoints, localTgtMat, 
reducedEpsilonAtReducedintegPoints, optimize = True) 


where 
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e reducediIntegrationWeights[g] = Ôg e R":, 

e reducedEpsilonAtReducedIntegPoints[1gj] = €& (w) (îs) € 
RIX, 

e localTgtMat[glm] =Ki.m (eú), Yn) (îe) e R"sx4xd, 


The object reducedEpsilonAtReducedIntegPoints isa third-order tensor 
containing the components of the strain tensor applied to the ROB and evaluated at 
the reduced quadrature points. This quantity is precomputed in the offline stage. The 
object LocalTgtMat is a third-order tensor containing the components of the local 
tangent matrix evaluated at the reduced quadrature points: this quantity is computed 
online by the constitutive law solver. 


4.2.4 Tutorials and Datasets 


Various tutorials detailing some physical uses cases and capabilities of genericROM 
are available online.? Namely, variants of the POD-ECM methods are presented for 
nonlinear structural mechanics and nonlinear transient thermal problems, featuring 
various boundary conditions. The datasets needed to run these tutorials are also 
provided. 

We detail here one tutorial,!° corresponding to a nonlinear structural mechanics 
case, as described in [1]. To run this tutorial, a Z-set!! license is required for the 
constitutive law solver (Z-mat). 


4.2.4.1 Description of the Physics Problem 


Consider a 1 m-wide cube, illustrated in Fig. 4.6. This cube is subjected to 


e a variable thermal loading, in the form of time- and space-dependent scalar field 
obtained from a previous thermal computation, see Fig. 4.6 (right), 

e acentrifugal effect generated by a rotation along the Z-axis, 

e atime- and space-dependent pressure filed applied on the X, face, 

e fixed displacement imposed along the X-axis on the Xo face, along the Y-axis on 
the Xo Yo edge and along the Z-axis on the XoY; Z; vertex. 


The modeling hypothesis are: quasistatic equilibrium equations for the deformable 
solid in small perturbations. The material is elastoviscoplastic with nonlinear hard- 
ening and a Norton flow. The constitutive law is specified by the .mat file. 


? https://genericrom.readthedocs.io/en/latest/_tutorials/index.html. 
10 https://genericrom.readthedocs.io/en/latest/_tutorials/Mechanical/MecaSequential/index.html. 
l http://www.zset-software.com. 
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Fig. 4.6 (left) Dual mesh of the test case with the accumulated plasticity filed at the end of the 
simulation, (right) mesh and maximal temperature loading field (genericROM documentation) 


4.2.4.2 Reduction Strategy 


This is a simple “reproducting use case”, where a reduced-order model is constructed 
to reproduce a high-fidelity solution and check its quality. This use case contains the 
modeling complexity of the high-pressure turbine blade cyclic extrapolation, except 
for the parallel computation in distributed memory and the fact that the turbine fea- 
tures two material (elastic and elastoviscoplastic). The number of simulation cycles 
can be seen as the input, while the outputs are the displacement field “U” and the 
accumulated plastic strain field p (named “evrcum’”). 

The reduced-order basis is generated by the snapshot-POD method “Compress- 
Data” from the “DataCompressors.FusedSnapshotPOD” module. 

The operator compression step is computed using the ECM (Empirical Cubature 
Method) “CompressOperator” from the “OperatorCompressors.Mechanical” mod- 
ule. 


4.2.4.3 Algorithm 


Snapshots are generated by Z-set and read from the Z-set format. The quality of 
the data compression is evaluated by computed the relative projection error of the 
high-fidelity solutions on the reduced-order basis: 


|e- E (e w E 
—____—_—— <1x 10”, (4.4) 


where Y“ denote the POD modes (or vectors of the reduced-order basis). The same 
quantity is considered from p. 

The quality of the complete procedure is evaluated by compute the relative £2- 
norm error between the high-fidelity solutions and the reduced ones: 
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N Uyyu 
ua YY; 


luli 


< 1x105, (4.5) 


with y;“ the coefficients of the reduced solution (or “reducedCoordinates”) computed 
by the reduced-order model. 


4.2.4.4 Code—Offline Stage 


List of imports required for the offline stage of this example: 


from genericROM.IO import ZsetMeshReader as ZMR 

from genericROM.IO import ZsetSolutionReader as ZSR 

from Mordicus.Containers import ProblemData as PD 

from Mordicus.Containers import CollectionProblemData as CPD 
from Mordicus.Containers import Solution as S 

from genericROM.FE import FETools as FT 

from genericROM.DataCompressors import FusedSnapshotPOD as SP 
from genericROM.OperatorCompressors import Mechanical 

from Mordicus.IO import StateIO as SIO 

import numpy as np 


Then, filename and dimensions related to the mesh and the solutions have to be 
declared, readers are initalized, and the mesh is read: 


folder = 
GetTestDataPath()+"Zset"+os.sep+"MecaSequential"+os.sep 


meshFileName = folder + "cube.geof" 

solutionFileName = folder + "cube.ut" 

meshReader = ZMR.ZsetMeshReader (meshFileName) 

solutionReader = ZSR.ZsetSolutionReader (solutionFileName) 

mesh = meshReader.ReadMesh () 

numberOfNodes = mesh.GetNumberOfNodes () 

numberOfIntegrationPoints = 
FT.ComputeNumberOfIntegrationPoints (mesh) 

nbeOfComponentsPrimal = 3 

nbeOfComponentsDual = 6 


Then, the part of the ECM algorithm depending only on the mesh is carried out: 


operatorPreCompressionData = 
Mechanical . PreCompressOperator (mesh) 


Then, the objects “Solution” are declared and populated with data from the pre- 
computed Z-set solutions: 


outputTimeSequence = solutionReader.ReadTimeSequenceFromSolutionFile() 
solutionU = §.Solution("U", nbeOfComponentsPrimal, numberOfNodes, primality 
= True) 
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solutionSigma = S.Solution("sigma", nbeOfComponentsDual, 
numberOfIntegrationPoints, primality = False) 
solutionEvrceum = S.Solution("evrcum", 1, numberOfIntegrationPoints, 
primality = False) 
for time in outputTimeSequence: 
solutionU.AddSnapshot (solutionReader.ReadSnapshot("U", time, 
nbeOfComponentsPrimal, primality=True), time) 
solutionSigma.AddSnapshot (solutionReader.ReadSnapshot ("sig", time, 
nbeOfComponentsDual, primality=False), time) 
solutionEvrcum.AddSnapshot (solutionReader.ReadSnapshotComponent ("evrcum", 
time, primality=False), time) 


Then, the objects “CollectionProblemData” and “ProblemData” are declared, 
which will enable to agregate the “Solution” objects previously constructed in a 
standard fashion in Mordicus: 


problemData = PD.ProblemData("MecaSequential") 

problemData.AddSolution(solutionv) 

problemData.AddSolution(solutionSigma) 

problemData.AddSolution(solutionEvrcum) 

collectionProblemData = CPD.CollectionProblemData () 

collectionProblemData.AddVariabilityAxis(’config’, str, 
description="dummy variability") 


collectionProblemData.DefineQuantity("U", "displacement", "m") 

collectionProblemData.DefineQuantity("sSigma", "stress", "Pa") 

collectionProblemData.DefineQuantity("evreum", "accumulated 
plasticity"; "t 


collectionProblemData.AddProblemData(problemData, 
config="case-1" 


Then, the L2(2) correlation operator between snapshots is computed (identified 
by “U”): 


snapshotCorrelationOperator = 
{"U":FT.ComputeL2ScalarProducMatrix(mesh, 3) } 


Then, using the snapshots-POD method, we compute the reduced-order basis for 
the solutions “U” with the L2 (Q) correlation operator, and for the solutions p without 
correlation operator (the default operator is the identity) as a precomputing step for 
the Gappy-POD reconstruction method on p. 


SP.CompressData(collectionProblemData, "U", 1.e-6, 
snapshotCorrelationOperator["U"]) 
SP.CompressData(collectionProblemData, "“evrcum", 1.e-6) 


Then, we compute the reduced coefficients (or “reducedCoordinates”’) by project- 
ing the high-fidelity snapshots onto the reduced-order basis: 


collectionProblemData.CompressSolutions("U", 
snapshotCorrelationOperator["U"]) 
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Notice that the two previous steps can be done in one by setting the attribute 
“compressSolutions = True” in the function “SP.CompressData”. Then, we verify 
the quality of the data compression on “U”: 


reducedOrderBasisU = collectionProblemData.GetReducedOrderBasis("U") 
CompressedSolutionU = solutionU.GetReducedCoordinates() 
compressionErrors = [] 


for t 


in outputTimeSequence: 


reconstructedCompressedSolution = np.dot(CompressedSolutionU[t], 


reducedOrderBasisuU) 


exactSolution = solutionU.GetSnapshot (t) 
norml2ExSol = np.linalg.norm(exactSolution) 


if 


norml2ExSol != 0: 

relError = 
np.linalg.norm(reconstructedCompressedSolution-exactSolution) 

/norml2ExSol 


else: 


relError = 
np.linalg.norm(reconstructedCompressedSolution-exactSolution) 
compressionErrors.append(relError) 


Then, 
scheme: 


we carry out the ECM algorithm to determine the reduced quadrature 


Mechanical .CompressOperator (collectionProblemData, 
operatorPreCompressionData, mesh, 1.e-5, 


listNameDualVarOutput = ["evrcum"], 


listNameDualVarGappyIndicesforECM = ["evrcum"]) 


Finally, at the end of the offline, the Modicus datamodel containing the results of 
this stage, is saved on disk in order to use it during the online stage. 


SIO.SaveState("collectionProblemData", collectionProblemData) 
SIO.SaveState("snapshotCorrelationOperator", 
snapshotCorrelationOperator) 


4.2.4.5 


Code—Online Stage 


List of imports required for the offline stage of this example: 


from 
from 
from 
from 
from 
from 
from 
from 
from 


genericROM.IO import ZsetInputReader as ZIR 
genericROM.IO import ZsetMeshReader as ZMR 

genericROM.IO import ZsetSolutionReader as ZSR 
genericROM.IO import ZsetSolutionWriter as ZSW 
Mordicus.Containers import ProblemData as PD 
Mordicus.Containers import Solution as S 

genericROM.FE import FETools as FT 
genericROM.OperatorCompressors import Mechanical as Meca 
Mordicus.IO import StateIO as SIO 


import numpy as np 
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First, data saved on disk at the end of the offline stage is read: 


collectionProblemData = SIO.LoadState("collectionProblemData" ) 
operatorCompressionDataMechanical = 
collectionProblemData.GetOperatorCompressionData("U") 
snapshotCorrelationOperator = 
SIO.LoadState("snapshotCorrelationOperator") 
reducedOrderBases = 
collectionProblemData.GetReducedOrderBases () 


Then, filename and dimensions related to the mesh and the solutions have to be 
declared and readers are initialized, in the same fashion as the offline stage. We 
mention here that the physical setting for the online stage has been taken identical 
to the one used in the offline stage (the folder “MecaSequential/”). In meaningful 
workflow, the physical setting for the online stage would be different. 


folder = 
GetTestDataPath()+"Zset"+os.sep+"MecaSequential"+os.sep 

inputFileName = folder + "cube.inp" 

inputReader = ZIR.ZsetInputReader (inputFileName) 

meshFileName = folder + "cube.geof" 


Then, the mesh is read (which is required when the variability is not parametrized): 


mesh = ZMR.ReadMesh (meshFileName) 


Then, an object “ProblemData” is defined, which will store the data computed 
during the online stage: 


onlineProblemData = PD.ProblemData("Online") 
onlineProblemData.SetDataFolder(os.path.relpath(folder, 
folderHandler.scriptFolder) ) 


Then, the temporal sequence and the constitutive law are read from the Z-Set 
input file. These are “fixed data” for the online resolution: 


timeSequence = inputReader.ReadInputTimeSequence() 

constitutiveLawsList = 
inputReader.ConstructConstitutiveLawsList () 

onlineProblemData.AddConstitutiveLaw(constitutiveLawsList) 


Then, the loadings and initial condition are read from the Z-Set input file and are 
reduced by projecting them onto the reduced-order basis: 


loadingList = inputReader.ConstructLoadingsList () 
onlineProblemData.AddLoading(loadingList) 
for loading in onlineProblemData.GetLoadingsForSolution("U"): 
loading.ReduceLoading(mesh, onlineProblemData, 
reducedOrderBases, operatorCompressionData) 
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initialCondition = inputReader.ConstructInitialCondition() 

onlineProblemData.SetInitialCondition(initialCondition) 

initialCondition.ReduceInitialSnapshot (reducedOrderBases, 
snapshotCorrelationOperator) 


Then, the reduced solution is computed in a nonintrusive fashion using a reduced 
Newton iterative algorithm for solving the reduced nonlinear system of equations at 
each time-step: 


onlineCompressedSolution = 
Meca.ComputeOnline(onlineProblemData, timeSequence, 
operatorCompressionDataMechanical, 1.e-8) 


Then, the reduced coefficients (or “reducedCoordinates”’) for the dual quantified 
of interest, here p, are computed using the online part of the Gappy-POD: 


onlineData = onlineProblemData.GetOnlineData("U") 
onlineEvrcumCompressedSolution, errorGappy = 
Meca.ReconstructDualQuantity("evrcum", 
operatorCompressionDataMechanical, 
onlineData, timeSequence = 
list (onlineCompressedSolution.keys() ) ) 


In order to compare the reduced solutions to the high-fidelity reference ones, 
“Solution” objects are created and populated with precomputed Z-set solutions: 


numberOfIntegrationPoints = 
FT.ComputeNumberOfIntegrationPoints (mesh) 
nbeOfComponentsPrimal = 3 
numberOfNodes = mesh.GetNumberOfNodes () 
solutionFileName = folder + "cube.ut" 
solutionReader = ZSR.ZsetSolutionReader (solutionFileName) 
outputTimeSequence = 
solutionReader .ReadTimeSequenceFromSolutionFile() 


solutionEvrcumExact = S.Solution("evrcum", 1, 
numberOfIntegrationPoints, primality = False) 
solutionUExact = S.Solution("U", nbeOfComponentsPrimal, 
numberOfNodes, primality = True) 
for t in outputTimeSequence: 
evrcum = solutionReader.ReadSnapshotComponent("evrcum", t, 
primality=False) 
solutionEvrcumExact.AddSnapshot (evreum, t) 
U = solutionReader.ReadSnapshot("U", t, 
nbeOfComponentsPrimal, primality=True) 
solutionUExact.AddSnapshot(U, t) 


“Solution” objects corresponding to the reduced solutions are constructed, popu- 
lated with the reduced coefficients (or “reducedCoordinates”’) computed by the online 
stage, and reconstructed on the complete mesh: 
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solutionEvrcumApprox = S.Solution("evrcum", 1, numberOfIntegrationPoints, 
primality = False) 

solutionEvrcumApprox. SetReducedCoordinates (onlineEvrcumCompressedSolution) 

solutionEvrcumApprox.UncompressSnapshots (reducedOrderBases ["evrcum" ] ) 

solutionUApprox = §S.Solution("U", nbeOfComponentsPrimal, numberOfNodes, 
primality = True) 

solutionUApprox. SetReducedCoordinates (onlineCompressedSolution) 

solutionUApprox.UncompressSnapshots (reducedOrderBases["U"] ) 


Then, we verify the quality of the reduced solutions “U” and “evrcum”: 


ROMErrorsU = [] 
ROMErrorsEvrcum = [] 
for t in outputTimeSequence: 
exactSolution = solutionEvrcumExact.GetSnapshotAtTime(t) 
approxSolution = solutionEvrcumApprox.GetSnapshotAtTime (t) 
norml2ExactSolution = np.linalg.norm(exactSolution) 
if norml2ExactSolution > 1.e-10: 
relError = 
np.linalg.norm(approxSolution-exactSolution) /norml2ExactSolution 
else: 
relError = np.linalg.norm(approxSolution-exactSolution) 
ROMErrorsEvrcum. append(relError) 


exactSolution = solutionUExact.GetSnapshotAtTime(t) 
approxSolution = solutionUApprox.GetSnapshotAtTime (t) 
norml2ExactSolution = np.linalg.norm(exactSolution) 
if norml2ExactSolution > 1.e-10: 

relError = 

np.linalg.norm(approxSolution-exactSolution) /norml2ExactSolution 

else: 

relError = np.linalg.norm(approxSolution-exactSolution) 
ROMErrorsU. append (relError) 


Finally, reduced predictions for “U” and “evrcum” are exported in the Z-set 
format: 


onlineProblemData.AddSolution(solutionUApprox) 

onlineProblemData.AddSolution(solutionEvrcumApprox) 

ZSW.WriteZsetSolution(mesh, meshFileName, "reduced", 
collectionProblemData, onlineProblemData, "U") 


4.2.4.6 Results 


A comparison between the reduced and reference high-fidelity solutions is illustrated 
in Fig. 4.7. 
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Fig. 4.7 Comparison between the reduced and reference high-fidelity solutions: (top) on the dis- 
placement “U”, (bottom) on the accumulated plasticity “p” (genericROM documentation) 
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Chapter 5 A) 
Industrial Application: Uncertainty get 
Quantification in Lifetime Prediction 

of Turbine Blades 


This chapter is a synthesis of the previous ones, since many introduced concepts are 
applied herein. The complete ROM-net workflow, described in Sect. 2.4.2 is applied 
to the quantification of the uncertainty of dual quantities (such as the accumulated 
plastic strain and the stress tensor) on an real-life turbine blade, generated by the 
uncertainty of the temperature loading field. The numerical experiments make use of 
the codes Mordicus and genericROM, introduced respectively in Sects. 4.1 and 4.2. 
The content of this chapter is inspired from our publication [9]. 

Computing the fatigue lifetime of one such blade requires simulating its behavior 
until the stabilization of the mechanical response, which last several weeks using 
Abaqus [23] because of the size of the mesh, the complexity of the constitutive 
equations, and the number of loading cycles in the transient regime. With such a 
computation time, uncertainty quantification with the Monte Carlo method is unaf- 
fordable. In addition, such simulations are too time-consuming to be integrated in 
design iterations, which limits them to the final verification steps, while the design 
process still relies on simplified models. Accelerating these complex simulations is 
a key challenge while maintaining a satisfying accuracy, as it would provide useful 
numerical tools to improve design processes and quantify the effect of the uncertain- 
ties on the environment of the system. 

Simulations are accelerated using a dictionary of reduced order models, with a 
classifier able to select which local reduced order model to be used for a new tem- 
perature loading. A dataset of 200 solutions is computed in a Finite Element approx- 
imation space of dimension in the order of the million, for various instances of the 
temperature field loading, in parallel in 7 days and 9 h on 48 cores. These solutions 
are computed over 11 time steps in the first cycle, using a scalable Adaptive Mul- 
tiPreconditioned FETI (AMPFETI) solver [3] in Z-set finite-element software [17]. 
The dataset is partitioned into two clusters using a k-medoids algorithm with a ROM- 
oriented dissimilarity measure in 5 min; the corresponding local ROMs, using POD 
data compression and ECM operator compression, are trained in 2 h and 30 min. An 
automatic reduced model recommendation procedure, allowing to decide which local 
ROM to be used for a new temperature loading, is trained in the form of a logistic 
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regression classifier in 16 min. A meta-model is used to reconstruct the dual quanti- 
ties of interest over the complete mesh from their values on the reduced integration 
points, in the form of a multi-task Lasso, which takes 1 h to train for 14 dual fields. 
The uncertainties on dual quantities of interest, such as the accumulated plastic strain 
and the stress tensor, are quantified by using our trained ROM-net on 1008 Monte 
Carlo draws of the temperature loading field in 2 h and 48 min, which corresponds to 
a speedup greater than 600 with respect to our highly optimized domain decomposi- 
tion AMPFETI solver. Expected values for the Von Mises stress and the accumulated 
plastic strain have 0.99-confidence intervals width of respectively 1.66% and 2.84% 
(relative to the corresponding prediction for the expected value). As a verification 
stage, 20 reference solutions are computed on new temperature loadings, and dual 
quantities of interest are predicted with relative accuracy in the order of 1% to 2%, 
while the location of the maximum value is perfectly predicted. 

In what follows, we describe the industrial dataset, the hypotheses of the model, 
and the objective of the present study. The proposed workflow for uncertainty quan- 
tification is then applied on this industrial configuration. 


5.1 Industrial Context 


The industrial test case of interest consists in predicting the mechanical behavior of 
a high-pressure (HP) turbine blade in an aircraft engine with uncertainties on the 
thermal loading. The industrial context and the models for the mechanical behavior 
and the thermal loading are presented, with a particular emphasis on the assumptions 
that have been made. For confidentiality reasons, mesh sizes and numerical values 
corresponding to the industrial dataset are not given, and the provided figures and 
plots do not contain any color map or physical numerical value. The accuracy of the 
predictions made by our methodology are given in the form of relative errors. 


5.1.1 Thermomechanical Fatigue of High-Pressure Turbine 
Blades 


High-pressure turbine blades are critical parts in an aircraft engine. Located down- 
stream of the combustion chamber, they are subjected to extreme thermomechanical 
loadings resulting from the combination of centrifugal forces, pressure loads, and hot 
turbulent fluid flows whose temperatures are higher than the material’s melting point. 
The repeating thermomechanical loading over time progressively damages the blades 
and leads to crack initiation under thermomechanical fatigue. Predicting the fatigue 
lifetime is crucial not only for safety reasons, but also for ecological issues, since 
reducing fuel consumption and improving the engine’s efficiency requires increasing 
the temperature of the gases leaving the combustion chamber. 
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High-pressure turbine blades are made of monocrystalline nickel-based super- 
alloys that have good mechanical properties at high temperatures. To reduce the 
temperature inside this material, the blades contain cooling channels where flows 
relatively fresh air coming from the compressor. In addition, the blade’s outer sur- 
face is protected by a thin thermal barrier coating. In spite of these advanced cooling 
technologies, the rotor blades undergo centrifugal forces at high temperatures, caus- 
ing inelastic strains. Under this cyclic thermomechanical loading repeated over the 
flights, the structure has a viscoplastic behavior and reaches a viscoplastic stabilized 
response, where the dissipated energy per cycle still has a nonzero value. This is called 
plastic shakedown, and leads to low-cycle fatigue. At cruise flight, the persistent cen- 
trifugal force applied at high temperature induces progressive (or time-dependent) 
inelastic deformations: this phenomenon is called creep. In addition, the difference 
between gas pressures on the extrados and the intrados of the blade generates bending 
effects. Environmental factors may also locally modify the chemical composition of 
the material, leading to its oxidation. As oxidized parts are more brittle, they facilitate 
crack initiation and growth. Thermal fatigue resulting from temperature gradients is 
another life-limiting factor. Temperature gradients make colder parts of the structure 
prevent the thermal expansion of hotter parts, creating thermal stresses. Due to their 
higher temperatures, the hot parts are more viscous and have a lower yield stress, 
which make them prone to develop inelastic strains in compression. When the tem- 
perature cools down after landing, tensile residual stresses appear in parts which 
were compressed at high temperatures and favor crack nucleation. Given the com- 
plex temperature field resulting from the internal cooling channels and the turbulent 
gas flow, thermal fatigue has a strong influence on the turbine blade’s lifetime. In 
particular, during transient regimes such as take-off, an important temperature gra- 
dient appears between the leading edge and the trailing edge of the blade, since the 
latter has a low thermal inertia due to its small thickness and thus warms up faster. 

In short, the behavior of a high pressure turbine blade results from a complex 
interaction between low-cycle fatigue, thermal fatigue, creep, and oxidation. Due to 
the cost and the complexity of experiments on parts of an aircraft engine, numerical 
simulations play a major role in the design of high-pressure turbine blades and 
their fatigue lifetime assessment. All this knowledge have been learned by scientist 
and engineers during last decades. In the proposed approach to machine learning for 
model order reduction, all this knowledge is preserved in local ROMs. It is even more 
than that, the uncertainty propagation comes to complete this valuable traditional 
knowledge. We do not expect from artificial intelligence to learn everything in our 
modeling process. 


5.1.2 Industrial Dataset and Objectives 


Figure 5.1 gives the geometry and the finite-element mesh of a real high-pressure 
turbine blade. The mesh is made of quadratic tetrahedral elements, and contains 
a number of nodes in the order of the million. The elasto-viscoplastic mechani- 
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Fig. 5.1 High-pressure turbine blade geometry and mesh (micro-perforations are not modeled) [9] 


cal behavior is described by a crystal plasticity model presented in the appendix 
of [9]. As explained above, Monte Carlo simulations using a commercial software 
as Abaqus are unaffordable. With the help of domain decomposition methods, the 
computation time can be reduced by solving equilibrium equations in parallel on 
different subdomains of the geometry. Using the implementation of the Adaptive 
MultiPreconditioned FETI solver [3] in Z-set finite-element software [17], the simu- 
lation of one single loading cycle of the HP turbine blade with 48 subdomains takes 
approximately 53 min. 

The objective is to use a ROM-net to quantify uncertainties on the mechanical 
behavior of the high-pressure turbine blade, given uncertainties on the thermal load- 
ing. The reduction of the computation time should enable Monte Carlo simulations 
for uncertainty quantification. In particular, we are not interested in predicting the 
state of the structure after a large number of flight-representative loading cycles. 
Only one cycle is simulated. Cyclic extrapolation of the behavior of a high-pressure 
turbine blade has been studied in [4] and is out of the scope of this section. 


5.1.3 Modeling Assumptions 


It is assumed that the heat produced or dissipated by mechanical phenomena has 
negligible effects in comparison with thermal conduction, which enables avoiding 


5.1 Industrial Context 75 


Fig. 5.2 Function w(t) w(t) 
defining one cycle for the 
rotation speed [9] 


strongly coupled thermomechanical simulations and running thermal and mechan- 
ical simulations separately instead. Under a weak thermomechanical coupling, the 
first step consists in solving the heat equation to determine the temperature field and 
its evolution over time. The temperature field history defines the thermal loading 
and is used to compute thermal strains and temperature-dependent material param- 
eters for the mechanical constitutive laws. Once the thermal loading is known, the 
temperature-dependent mechanical problem must be solved in order to predict the 
mechanical response of the structure. 

The thermomechanical loading applied to the high-pressure turbine blade during 
its whole life is modeled as a cyclic loading, with one cycle being equivalent to one 
flight. The rotation speed of the turbine’s rotor is proportional to a periodic function 
of time w(t) whose evolution over one period (or cycle, see Fig. 5.2) is representative 
of one flight with its three main regimes, namely take-off, cruise, and landing. The 
period (or duration of one cycle) is denoted by t+. The rotation speed between flights 
k and k + 1 is zero, which means that w(kt,) = 0 for any integer k. The rotation 
speed w(t) is scaled so that its maximum is 1. 

Let Q C R? denote the solid body representing the high-pressure turbine blade, 
with dQ denoting its outer surface. Let QP C dQ be the surface corresponding to 
the intrados and extrados. The thermal loading is defined as: 


VEEQ, WeR,, TE,t)=(— w(t))To + w(t)Tma (f), (3.1) 
where Tọ = 293 K and Tmax 1s the temperature field obtained when the rotation speed 
reaches its maximum. This field Tmax is obtained either by an aerothermal simulation 


or by a stochastic model, as explained later. Similarly, the pressure load applied on 
dQ? reads: 


VE cI, WER, pE, t= (1 — a(t) po? + a(t) p?2(&), 6D 
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where pł? = 1 atm is the atmospheric pressure at sea level, and where p?®, is the 
pressure field obtained when the rotation speed reaches its maximum. The clamping 
of the blade’s fir-tree foot on the rotor disk is modeled by displacements boundary 
conditions that are not detailed here. 

Small geometric details of the structure have been removed to simplify the geome- 
try. Nonetheless, the main cooling channels are considered. The effects of the thermal 
barrier coating (TBC) have been integrated in aerothermal simulations, but the TBC 
is not considered in the mechanical simulation although its damage locally increases 
the temperature in the nickel-based superalloy and thus affects the fatigue resistance 
of the structure. Additional centrifugal effects due to the TBC are not taken into 
account. 

The predicted mechanical response of the structure depends on many different 
factors. Below is a nonexhaustive list of influential factors that are possible sources 
of uncertainties in the numerical simulation: 


e Thermal loading: The viscoplastic behavior of the nickel-based superalloy is 
very sensitive to the temperature field and its gradients. However, the temperature 
field is not accurately known because of the impossibility of validating numerical 
predictions experimentally. Indeed, temperature-sensitive paints are accurate to 
within 50 K only, and they do not capture a real surface temperature field since 
they measure the maximum temperature reached locally during the experiment. 

e Crystal orientation: Because of the complexity of the manufacturing process of 
monocrystalline blades, the orientation of the crystal is not perfectly controlled. As 
the superalloy has anisotropic mechanical properties, defaults in crystal orientation 
highly affect the location of damaged zones in the structure. 

e Mechanical loading: The centrifugal forces are well known because they are 
related to the rotation speed that is easy to measure. On the contrary, pressure loads 
are uncertain because of the turbulent nature of the incoming fluid flow. However, 
the effects of pressure loads uncertainties on the mechanical response are less 
significant than those of the thermal loading and crystal orientation uncertainties. 

e Constitutive laws: Uncertainties on the choice of the constitutive model, the 
relevance of the modeling assumptions, and the values of the calibrated parameters 
involved in the constitutive equations also influence the results of the numerical 
simulations. 


For simplification purposes, the only source of uncertainty that is considered 
in this work is the thermal loading. The equations of the mechanical problem are 
then seen as parametrized equations, where the parameter is the temperature field 
Tmax (See Eq. (5.1)) obtained when the rotation speed reaches its maximum value. 
The dimension of the parameter space is then the number of nodes in the finite- 
element mesh. The mechanical loading is assumed to be deterministic. With the 
crystal orientation, the constitutive laws and their parameters (or coefficients), they 
are considered as known data describing the context of the study and given by experts. 
For details on the constitutive law model, we refer to the appendix of [9]. 
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5.1.4 Stochastic Model for the Thermal Loading 


A stochastic model is required to take into account the uncertainties on the thermal 
loading. Given the definition of the thermal loading in Eq. (5.1), we only need to 
model uncertainties in space through the field Tmax obtained when the rotation speed 
reaches its maximum value. The random temperature fields must satisfy some con- 
straints: they must satisfy the heat equation, and they must not take values out of the 
interval [0 K; Tet], where Tmeit is the melting point of the superalloy. These random 
fields are obtained by adding random fluctuations to a reference temperature field, 
see Fig. 5.3. The reference field and comes from aerothermal simulations run with 
the software Ansys Fluent.' The data-generating distribution is defined as a Gaussian 
mixture model made of two Gaussian distributions with the same covariance func- 
tion but with distinct means, and with a prior probability of 0.5 for each Gaussian 
distribution. The Gaussian distributions are obtained by taking the four first eigen- 
functions of the covariance function (see Karhunen-Loéve expansion [16]), with a 
standard deviation of 15 K. Therefore, realizations of the random temperature field 
read: 


4 
VE EQ, T(E) = Telé) + Yo STE) + D> Ti 8T), (5.3) 


i=1 


where Tef is the reference field, 57> is a temperature perturbation at the trailing 
edge whose maximum value is 50 K, {ôT; }<i<4 are fluctuation modes, Yo is a ran- 
dom variable following the Bernoulli distribution with parameter 0.5, and {Y;}1<;<4 
are independent and identically distributed random variables following the standard 
normal distribution N(0, 1). The variable Yo is also independent of the other vari- 
ables Y;. The different fields involved in Eq. (5.3) can be visualized in Fig. 5.3. 
Equation (5.3) defines a mixture distribution with two Gaussian distributions whose 
means are Tyrer and Tief + 67. We voluntarily define this mixture distribution with 
ôT adding 50 K in a critical zone of the turbine blade in order to check that our clus- 
ter analysis can successfully detect two relevant clusters, i.e., one for fields obtained 
with Yo(@) = 0 and one for fields obtained with Yo(@) = 1. Indeed, the temperature 
perturbation 57 is expected to significantly modify the mechanical response of the 
high-pressure turbine blade. All the fields {57;}o<j<4 satisfy the steady heat equa- 
tion like Tef, which ensures that the random fields always satisfy the heat equation 
under the assumption of a linear thermal behavior. For nonlinear thermal behaviors, 
Eq. (5.3) would define surface temperature fields that would be used as Dirichlet 
boundary conditions for the computation of bulk temperature fields. The assumption 
of a linear thermal behavior is adopted here to avoid solving the heat equation for 
every realization of the random temperature field. 

Let us now give more details about the construction of the fluctuation modes 
{ôT; }<i<4. First, surface fluctuation modes are computed on the boundary 0Q using 
the method given in [21] for the construction of random fields on a curved surface. 


l https://www.ansys.com/products/fluids/ansys-fluent. 
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Fig. 5.3 Reference temperature field (on the left), temperature perturbation at the trailing edge 


(field 0 = 579), and fluctuation modes (fields 1 to 4). The fluctuations in the fourth mode are located 
inside the blade, in the cooling channels [9] 


The correlation function is defined as a function of the geodesic distance dg along 


the surface dQ: ; 
dg(&, £) ) 


d? (5.4) 


pÈ, é’) = exp (- 


where de. is acorrelation length. Geodesic distances are computed using the algorithm 
described in [18, 22] and implemented in the Python library gdist.2 A covariance 
matrix is built by evaluating the correlation function on pairs of nodes of the outer 
surface of the finite-element mesh, and multiplying the correlation by the constant 
variance. The four surface modes are then obtained by finding the four eigenvectors 
corresponding to the largest eigenvalues of the covariance matrix. The steady heat 
equation with Dirichlet boundary conditions is solved for each of these surface modes 
to derive the 3D fluctuation modes, using Z-set [17] finite-element solver. The Python 
library BasicTools* developed by SafranTech is used to read the finite-element mesh 
and write the temperature fields in a format that can be used for simulations on Z-set. 


? https://pypi.org/project/gdist/. 
3 https://gitlab.com/drti/basic-tools. 
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5.2 ROM-net Based Uncertainty Quantification Applied 
to an Industrial High-Pressure Turbine Blade 


This section develops the different stages of the ROM-net for the industrial test case 
presented in the previous section. Given our budget of 200 high-fidelity simula- 
tions, a dictionary containing two local ROMs is constructed using our clustering 
procedure. A logistic regression classifier is trained for automatic model recommen- 
dation using information identified by feature selection, followed by an alternative 
to the Gappy-POD for full-field reconstruction of dual quantities. Then, the results 
of the uncertainty quantification procedure are presented. Finally the accuracy of the 
ROM.-net is validated using simulations for new temperature loadings. 


5.2.1 Design of Numerical Experiments 


Given the computational cost of high-fidelity mechanical simulations of the high- 
pressure turbine blade, the training data are sampled from the stochastic model for 
the thermal loading using a design of experiments. Our computational budget corre- 
sponds to 200 high-fidelity simulations, so a database of 200 temperature fields must 
be built. This database includes two separate datasets coming from two independent 
DoEs: 


The first dataset is built from a Maximum Projection LHS design (MaxProj LHS 
DoE [13]) and contains 80 points. This dataset will be used for the construction 
of the dictionary of local ROMs via clustering. The MaxProj LHS DoE has good 
space-filling properties on projections onto subspaces of any dimension. 

The second dataset is built from a Sobol’ sequence (Sobol’ DoE) of 120 points. 
Using a suboptimal DoE method ensures that this second dataset is different and 
independent from the first one. The lower quality of this dataset with respect to 
the first one is compensated by its larger population. This dataset will be used for 
learning tasks requiring more training examples than the construction of the local 
ROMs, namely the classification task for automatic model recommendation, and 
the training of cluster-specific surrogate models for the reconstruction of full fields 
from hyper-reduced predictions on a reduced-integration domain. These surrogate 
models (Gappy surrogates) replace the Gappy-POD [11] method that is commonly 
used in hyper-reduced simulations to retrieve dual variables on the whole mesh. 


These DoEs are built with the platform Lagun.* The fact that these two datasets 
come from two separate DoEs is beneficial: as each of them is supposed to have 
good space-filling properties, they are both representative of the possible thermal 
loading and can therefore be used to define a training set and a test set for a given 
learning task. For instance, the classifier trained on the Sobol’ DoE can be tested 


4 https://gitlab.com/drti/lagun. 
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Fig. 5.4 Visualization of the MaxProj LHS DoE. The marginal distributions are represented on the 
diagonal. The 5D DoE is projected on 2D subspaces for visualization purposes, in order to check 
space-filling properties in 2D [9] 


on the MaxProj LHS DoE. The local ROMs built from snapshots belonging to the 
MaxProj LHS DoE can make predictions on the Sobol’ DoE that will be used for the 
training of the Gappy surrogates, which is relevant since the Gappy surrogates are 
supposed to analyze ROM predictions on new unseen data in the exploitation phase. 

Drawing random temperature fields as defined in Eq. (5.3) requires sampling data 
from the random variables {Y;}o<;<4, where Yo follows the Bernoulli distribution 
with parameter 0.5 and the variables Y; for i € [[1; 4] are independent standard 
normal variables and independent of Yo. Both DoE methods (Maximum Projection 
LHS and Sobol’ sequence) generate point clouds with a uniform distribution in 
the unit hypercube. Figures 5.4 and 5.5 show the projections onto 2-dimensional 
subspaces of the 5D point clouds used to build our datasets. The marginal distributions 
are plotted to check that they well approximate the uniform distribution. These point 
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Fig. 5.5 Visualization of the Sobol’ DoE. The marginal distributions are represented on the diag- 
onal. The 5D DoE is projected on 2D subspaces for visualization purposes, in order to check 
space-filling properties in 2D [9] 


clouds, considered as samples of a random vector (Xo, X1, X2, X3, X4) following the 
uniform distribution on the unit hypercube, are transformed into realizations of the 
random vector (Yo, Y1, Y2, Y3, Y4) using the following transformations: 


Yo = 1%>1⁄2 and Wie fl 4], =F"), (5.5) 


where F~! is the inverse of the cumulative distribution function of the standard 
normal distribution. The resulting samples define the MaxProj dataset and the Sobol’ 
dataset of random temperature fields, using Eq. (5.3). Each temperature field defines 
a thermal loading, using Eq. (5.1). The 200 corresponding mechanical problems are 
solved for one loading cycle with the finite-element software Z-set [17] with the 
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domain decomposition method described in [3], with 48 subdomains. The average 
computation time for one simulation is 53 min. 


5.2.2 ROM Dictionary Construction 


The 80 simulations associated to the MaxProj dataset are used as clustering data. 
Loading all the simulation data and computing the pairwise ROM-oriented dissim- 
ilarities takes about 5 min. The ROM-oriented dissimilarity defined in [7, Defini- 
tion 3.11] is computed with n = 1, i.e., each simulation is represented by one field. 
The dataset is partitioned into two clusters using our implementation of PAM [14, 15] 
k-medoids algorithm, with 10 different random initializations for the medoids. The 
clustering results can be visualized using Multidimensional Scaling (MDS) [2]. MDS 
is an information visualization method which consists in finding a low-dimensional 
dataset Zo whose matrix of Euclidean distances d (Zo) is an approximation of the true 
dissimilarity matrix ô. To that end, a cost function called stress function is minimized 
with respect to Z: 


Zo = aoe (¢(Z; 6)) = ag min Or — di; (Z)? : (5.6) 


i<j 


This minimization problem is solved with the algorithm Scaling by MAjoriz- 
ing a COmplicated Function (SMACOF, [10]) implemented in Scikit-Learn [19]. 
Figure 5.6 show the clusters on the MDS representations with the goal-oriented vari- 
ants of the ROM-oriented dissimilarity measure, applied on the accumulated plastic 
strain pim. The clustering results is compared with the expected clusters corre- 
sponding to Yo = 0 and Yo = 1, the latter corresponds to the perturbation ôTo being 
activated. The obtained clusters almost correspond to the expected ones, with only 4 
points with wrong labels out of 80, which quantifies the ability of the ROM-oriented 
dissimilarity measure on the accumulated plasticity to infer the correct value of Yo. 
The medoids of the two clusters are given in Fig. 5.7. Cluster 0 contains temperature 
fields for which Yo = 1, while cluster 1 contains fields for which Yo = 0. It can be 
observed that the quantity of interest clearly differs from one cluster to the other, 
while the differences are hardly visible on the displacement field. The displacement 
field combines deformations associated to different phenomena (thermal expansion, 
elastic strains, viscoplastic strains) that are not necessarily related to damage in the 
structure, which could explain why the quantity of interest p?,,, seems to be more 
appropriate for clustering in this example. 

The simulations used for the clustering procedure can directly provide snapshots 
for the construction of the local ROMs. To control the duration of their training, 
only 20 simulations are selected to provide snapshots for the each local ROMs, 
which represents half of the clusters’ populations. These simulations are selected in 
a maximin greedy approach starting from the medoid (see [8, Algorithm 2, Stage 2] 
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Fig. 5.6 MDS representation of the clustering results using the ROM-oriented dissimilarity measure 
on the quantity of interest pfam (goal-oriented variant). On the left, the colors correspond to the 
expected clusters. On the right, the colors correspond to the clusters identified by the clustering 
algorithm. The positions of the labels 0 and 1 coincide with the positions of the clusters’ medoids. 
The MDS relative error ¢ (Zo; 5)/¢ (0; ô) is 12% [9] 


Fig. 5.7 The 3 fields on the left correspond to the medoid of cluster 0, and those on the right 
correspond to the medoid of cluster 1. The fields in the first and the third columns show the 
differences between the medoids’ temperature fields and the reference temperature field Tref (the 
scale is truncated for the first field). The second and the fourth columns show the displacement 
magnitude field ./u.u (top) and the quantity of interest p? m (bottom) [9] 
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Fig. 5.8 MDS 
representation of the 
clustering results. Orange 
points represent the 
snapshots selected for cluster 
0, while the light blue points 
represent the snapshots 
selected for cluster 1. For 
each cluster, the snapshots 
are selected by a maximin 
procedure starting from the 
medoid [9] 


for a example of maximin selection). Figure 5.8 shows which simulations have been 
selected for the construction of the local ROMs. 

The local ROMs are built following the methodology described in Sect. 2.3. The 
snapshot-POD and the ECM are done in parallel with shared memory on 24 cores. The 
tolerance for the snapshot-POD is set to 1078 for the displacement field, and to 1074 
for dual variables (the quantity of interest p? m and the six components of the stress 
tensor). The POD bases for the dual variables will be used for their reconstruction 
with the Gappy surrogates. The tolerance for the ECM is set to 5 x 1074. Locals 
ROMs each contain 18 displacement modes and between 8 and 13 modes for stress 
components. The first and second local ROMs contain respectively 10 and 12 modes 
for the quantity of interest p? m- The ECM selects 506 (resp. 510) integration points 
for the reduced-integration domain of ROM 0 (resp. 1). Building one local ROM 
takes approximately 2 h and 30 min. 


5.2.3 Automatic Model Recommendation 


In this section, a classifier is trained for the automatic model recommendation task. 
The 120 temperature fields coming from the Sobol’ dataset are used as training data 
for the classifier. Their labels are determined by finding their closest medoid in terms 
of the ROM-oriented dissimilarity measure. Hence, for each temperature field of the 
Sobol’ dataset, two dissimilarities are computed: one with the medoid of the first 
cluster, and one with the medoid of the second cluster. Once trained, the classifier 
can be evaluated on the 80 labelled temperature fields of the MaxProj dataset. 

Each temperature field is discretized on the finite-element mesh, which contains 
in the order of the million nodes. To reduce the dimension of the input space and facil- 
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Fig. 5.9 Feature selection 
results. The kriging 
metamodel for redundancy 
terms is represented by the 
red curve and built from 800 
true redundancy terms (blue 
points). The elements 
containing the selected nodes 
are represented in the turbine 
blade geometry [9] 


Mutual Information 


Distance 


itate the training phase of the classifier, we apply the geostatistical mMRMR feature 
selection algorithm described in [8, Algorithm 1] on data from the Sobol’ dataset. 
First, 800 pairs of nodes are selected in the mesh, which takes 18 s. The 800 cor- 
responding redundancy terms are computed with Scikit-Learn [19] in less than 3 s. 
Figure 5.9 plots the values of these redundancy terms versus the Euclidean distance 
between the nodes. We observe that the correlation between the redundancy mutual 
information terms and the distance between the nodes is poor, with a lot of noise. 
This can be due to the fact that the random temperature fields have been built using 
Gaussian random fields on the outer surface with an isotropic correlation function 
depending on the geodesic distance along the surface rather than the Euclidean dis- 
tance. Since the turbine blade is a relatively thin structure, two nodes, one on the 
intrados and another one on the extrados, can be close to each other in the Euclidean 
distance, but with totally uncorrelated temperature fluctuations because of the large 
geodesic distance separating them. On the contrary, two points on the same side of 
the turbine blade can have correlated temperature variations while being separated 
by a Euclidean distance in the order of the blade’s thickness. The length of the mutual 
information’s high-variance regime seems to correspond to the blade’s chord, which 
supports this explanation. The thinness of the turbine blade induces anisotropy in the 
correlation function of the bulk Gaussian random field defining the thermal loading, 
which implies an anisotropic behavior of the mutual information according to [8, 
Property 1]. The use of a local temperature perturbation 57> in conjunction with 
fluctuation modes having larger length scales may also partially explain the large 
variance of redundancy terms. Nonetheless, it remains clear that redundancy terms 
are smaller as for large distances. This trend is captured by a kriging metamodel 
(Gaussian process regression) trained with Scikit-Learn in a few seconds, with a 
sum-kernel involving the Matérn kernel with parameter 5/2 (to get a continuous and 
twice differentiable metamodel) and length scale 1, and a white kernel to estimate 
the noise level of the signal. The curve of the metamodel is given in Fig. 5.9. Then, 
for each node of the finite-element mesh, the mutual information with the label vari- 
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Table 5.1 Classification results 


Class Precision Recall Fl-score Support 
0 0.9744 1.0000 0.9870 38 
1 1.0000 0.9762 0.9880 42 
Accuracy - - 0.9875 80 
Macro avg 0.9872 0.9881 0.9875 80 
Weighted avg 0.9878 0.9875 0.9875 80 


able is computed. The computations of these relevance terms (in the order of the 
million terms) are distributed between 280 cores, which gives a total computation 
time of 15 min. Among these features, 5, 986 features are preselected by discarding 
those with a relevance mutual information lower than 0.05. The geostatistical mRMR 
selects 11 features in 42 s. The corresponding nodes in the finite-element mesh can 
be visualized in Fig. 5.9. 


Remark 5.1 The metamodel for redundancy terms could be improved by defining 
it as a function of the precomputed geodesic distances along the outer surface rather 
than the Euclidean distances. Each finite-element node would be associated to its 
nearest neighbor on the outer surface before computing the approximate mutual 
information from geodesic distances. 


The classifier is trained on the Sobol’ dataset, using the values of the temperature 
fields at the 11 nodes identified by the feature selection algorithm. The classifier is 
a logistic regression [1, 5, 6] with elastic net regularization [24] implemented in 
Scikit-Learn. The two hyperparameters involved in the elastic net regularization are 
calibrated using 5-fold cross-validation, giving a value of 0.001 for the inverse of the 
regularization strength, and 0.4 for the weight of the L! penalty term (and thus 0.6 for 
the L? penalty term). Due to the L! penalty term, the classifier only uses 5 features 
among the 11 input features. The classifier’s accuracy, evaluated on the MaxProj 
dataset to use new unseen data, reaches 98.75%. The confusion matrix indicates that 
100% of the test examples belonging to class 0 have been correctly labeled, and that 
2.38% of the test examples belonging to class 1 have been misclassified. Table 5.1 
summarizes the values of precision, recall and Fl-score on test data. 


5.2.4 Surrogate Model for Gappy Reconstruction 


When using hyper-reduction, the ROM calls the constitutive equations solver only 
at the integration points belonging to the reduced-integration domain. It is recalled 
that the ECM selected 506 (resp. 510) integration points for the reduced-integration 
domain of ROM 0 (resp. 1), and that the finite-element mesh initially contains a 
number of integration points in the order of the million. Therefore, after a reduced 
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simulation, dual variables defined at integration points are known only at integra- 
tion points of the reduced-integration domain. To retrieve the full field, the Gappy- 
POD [11] finds the coefficients in the POD basis that minimize the squared error 
between the reconstructed field and the ROM predictions on the reduced-integration 
domain. This minimization problem defines the POD coefficients as a linear function 
of the predicted values on the reduced-integration domain. Although these coeffi- 
cients are optimal in the least squares sense, they can be biased by the errors made by 
the ROM. To alleviate this problem, we propose to replace the common Gappy-POD 
procedure by a metamodel or Gappy surrogate. The inputs and the outputs of the 
Gappy surrogate are the same as for the Gappy-POD: the input is a vector containing 
the values of a dual variable on the reduced-integration domain, and the output is a 
vector containing the optimal coefficients in the POD basis. One Gappy surrogate 
must be built for each dual variable of interest: in our case, 7 surrogate models per 
cluster are required, namely one for the quantity of interest p? m and one for every 
component of the Cauchy stress tensor. 

The training data for these Gappy surrogates are obtained by running reduced 
simulations with the local ROMs, using the thermal loadings of the Sobol’ dataset. 
Indeed, the two local ROMs have been built on the MaxProj dataset, therefore thermal 
loadings of the Sobol’ dataset can play the role of test data for the ROMs. For each 
thermal loading in the Sobol’ dataset, the true high-fidelity solution is already known 
since it has been computed to provide training data for the classifier. In addition, 
the exact labels for these thermal loadings are known, which means that we know 
which local ROM to choose for each thermal loading of the Sobol’ dataset. Given 
ROM predictions on the reduced-integration domain, the optimal coefficients in the 
POD basis are given by the projections of the true prediction made by the high- 
fidelity model (the finite-element model) onto the POD modes. This provides the 
true outputs for the Sobol’ dataset, which can then be used as a training set for the 
Gappy surrogates. 

Given the high-dimensionality of the input data (there are more than 500 integra- 
tion points in the reduced-integration domains) with respect to the number of training 
examples (120 examples), a multi-task Lasso metamodel is used. The hyperparameter 
controlling the regularization strength is optimized by 5-fold cross-validation. Train- 
ing the 14 Gappy surrogates (7 for each cluster) takes 1 h. The Gappy surrogates select 
between 8% and 18% of the integration points in the reduced-integration domains, 
due to the L! regularization. The mean cross-validated coefficients of determination 
are 0.9637 (resp. 0.8935) for the quantity of interest for cluster 0 (resp. cluster 1), and 
range from 0.9404 to 0.9938 for stress components. These satisfying results mean 
that it is not required to train a kriging metamodel with the variables selected by 
Lasso to get nonlinear Gappy surrogates. The Gappy surrogates are then linear, just 
as the Gappy-POD. 

The accuracy gains provided by the Gappy surrogates with respect to classical 
Gappy-POD on the present industrial case is investigated in Fig. 5.10. Here, 24 
high-fidelity simulation in the first cluster are computed (with Yo = 0) as reference, 
and Gappy surrogates (using two meta models Lasso and ElasticNet) and classical 
Gappy-POD are computed using ROM 0. For both variants of meta models and 
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Fig. 5.10 Mean over the complete mesh of dual quantities of interest: accumulated plastic strain 
Peum (left) and component 33 of the stress tensor 033 (right), plotted as points where the x-coordinate 
is the reference value, and the y-coordinate is the considered reduced prediction [9] 


both quantities of interest: accumulated plastic strain and the component 33 of the 


stress tensor, Gappy surrogates provides more accurate predictions than the classical 
Gappy-POD. 


Remark 5.2 In this strategy, the local ROMs solve the equations of the mechanical 
problem, which enables using linear surrogate models to reconstruct dual variables. 
Using surrogate models from scratch instead of local ROMs would have been more 
difficult, given the nonlinearities of this mechanical problem and the lack of training 
data for regression. In addition, such surrogate models would require a parametriza- 
tion of the input temperature fields, whereas the local ROMs use the exact values 
of the temperature fields on the RID without assuming any model for the thermal 
loading. 


The dictionary-based ROM-net used for mechanical simulations of the high- 
pressure turbine blade is made of a dictionary of two local hyper-reduced order 
models and a logistic regression classifier. The classifier analyzes the values of the 
input temperature field at 11 nodes only, identified by our feature selection strategy. 
For a given thermal loading in the exploitation phase, after the reduced simulation 
with the local ROM recommended by the classifier, linear cluster-specific Gappy sur- 
rogates reconstruct the full dual fields (quantity of interest and stress components) 
from their predicted values on the reduced-integration domain. 


5.2.5 Uncertainty Quantification Results 


Once trained, the ROM-net can be applied for the quantification of uncertainties on 
the mechanical behavior of the HP turbine blade resulting from the uncertainties 
on the thermal loading. Since the ROM-net online operations can be performed 
sequentially on one single core, 24 cores are used in order to compute the solution 
for 24 thermal loadings at once. This way, 42 batches of 24 Monte Carlo simulations 
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{ 


Fig. 5.11 Histograms and probability density functions of the quantities of interest p?,,,, (left) and 
Teq (right) [9] 


Table 5.2 Widths of the confidence intervals (CI) for the expectations, expressed as percentages 
of the estimated expectations 


Estimated variable Confidence level Relative CI width (%) 
elPeum 0.95 2.16 
elPeum! 0.99 2.84 
e€|Feq] 0.95 1.26 
elFeq] 0.99 1.66 


are run in 2 h and 48 min. The 1008 thermal loadings used for this study are generated 
by randomly sampling points from the uniform distribution on the 5D unit hypercube 
and applying the transformation given in Eq. (5.5). 

The expected values of Poum and Oeq are estimated with the empirical means 


Za = 1 ‘1 Zi, where Z; are the corresponding samples. The variances of p%,,,,, and 


= p i i 2 1 ” zs \2 
Oeq are computed using the unbiased sample variance $f = Raa J (Zi — Zn) 3 
n— 
i=1 


The Central Limit Theorem gives asymptotic confidence intervals for the expected 
values: for all œ €]0; 1[, 


I, = |Z» -osy S2/n; Zn + b1-¢y/ syn |, (5.7) 


where ¢, denotes the quantile of order r of the standard normal distribution N (0, 1), 
and I, is an asymptotic confidence interval with confidence level 1 — œ for the 
expectation n: lim, +40 P(N € In) = 1 — a. The widths of the confidence intervals 
are expressed as a percentage of the estimated value for the expectations in Table 5.2. 

The probability density functions of the quantities of interest can be estimated 
using Gaussian kernel density estimation (see Sect. 6.6.1. of [12]). Figure 5.11 gives 
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simulations with 24 cores) 


Fig. 5.12 Workflow for the ROM-net methology applied to the considered industrial setting [9] 


the histograms and estimated distributions for p®,,,, and Og. The shapes of these dis- 
tributions highly depend on the assumptions made for the stochastic thermal loading. 


5.2.6 Workflow 


Figure 5.12 provides an illustration of the workflow and the computational time of 
each step presented above. 


5.2.7 Verification 


For verification purposes, the accuracy of the ROM-net is evaluated on 20 Monte 
Carlo simulations with 20 new thermal loadings. These thermal loadings are gen- 
erated by randomly sampling points from the uniform distribution on the 5D unit 
hypercube, and applying the transformation given in Eq. (5.5). The reduced simula- 


5.2 ROM-net Based Uncertainty Quantification Applied to an Industrial ... 91 


tions are run on single cores. The total computation time for generating a new thermal 
loading on the fly, selecting the corresponding reduced model, running one reduced 
simulation and reconstructing the quantities of interest is 4 min on average. As a 
comparison, one single high-fidelity simulation with Z-set [17] with 48 subdomains 
takes 53 min, which implies that the ROM-net computes 13.25 times faster. However, 
one high-fidelity simulation requires 48 cores for domain decomposition, whereas 
the ROM-net works on one single core. Hence, using 48 cores to run 48 reduced 
simulations in parallel, 636 reduced simulations can be computed in 53 min with the 
ROM-net, while the high-fidelity model only runs one simulation. In addition to the 
acceleration of numerical simulations, energy consumption is reduced by a factor of 
636 in the exploitation phase. In spite of the fast development of high-performance 
computing, numerical methods computing approximate solutions at reduced compu- 
tational resources and time are particularly important for many-query problems such 
as uncertainty quantification, where the intensive use of computational resources 
is a major concern. Model order reduction and ROM-nets play a prominent role 
toward green numerical simulations [20]. Of course, the number of simulations in 
the exploitation phase must be large enough to compensate the efforts made in the 
training phase, like in any machine learning or model order reduction problem. 

Figures 5.13 and 5.14 show the results for two simulations belonging to cluster 
O and cluster 1 respectively. These figures give the difference between the current 
temperature field as the reference one, i.e., the field T — Tref, and the resulting vari- 
ations of the quantity of interest predicted by the ROM-net and the high-fidelity 
model, i.e., po? ROM(T) — p?HF(T..¢) and p?HF(T) — p?-HF(T,.¢). The signs and the 
positions of the variations of the quantity of interest seem to be quite well predicted 
by the ROM-net. 

Let us introduce a zone of interest Q defined by all of the integration points 
at which p,m is higher than 0.4 x max p%,,(&) for the thermal loading defined by 
Tet + 679. This zone of interest contains 209 integration points. The values of the 
variables p? m and eq averaged over Q’ are denoted by Plum and Gq. Table 5.3 gives 
different indicators quantifying the errors made by the ROM-net: the L? relative 
errors on the whole domain Q and on the zone of interest Q’, the L® relative errors 
on Q and Q’, the relative errors on Peum and Geq, and the errors on the locations of 
the points where the fields p@,,,, and oeg reach their maxima. All the relative errors 
remain in the order of 1% or 2%, which validates the methodology. In addition, the 
ROM-net perfectly predicts the position of the critical points at which pfam and eq 
reach their maxima. Figure 5.15 shows errors on the quantities of interest. 
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Fig. 5.13 Comparison between high-fidelity predictions (middle column) and ROM-net’s pre- 
dictions (right-hand column). The field on the left represents the difference between the current 
temperature field (belonging to cluster 0) and the reference one. The other fields correspond to 
the increments of the quantity of interest pf m with respect to its reference state obtained with the 


reference temperature field [9] 


Table 5.3 Error indicators for the evaluation of the ROM-net on 20 new thermal loadings 


Error indicator 


Errors on pein (%) 


Errors on eq (%) 


Mean L? relative error on Q 1.14 0.84 
Mean L? relative error on Q’ | 0.75 1.46 
Mean L® relative erroron Q | 1.11 1.09 
Mean L% relative error on Q’ | 1.05 2.60 
Mean rel. err. on value 0.50 0.89 
averaged over Q’ 

Mean distance between 0 0 


maxima 
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Fig. 5.14 Comparison between high-fidelity predictions (middle column) and ROM-net’s pre- 
dictions (right-hand column). The field on the left represents the difference between the current 
temperature field (belonging to cluster 1) and the reference one. The other fields correspond to 
the increments of the quantity of interest p?,,, with respect to its reference state obtained with the 
reference temperature field [9] 


Fig. 5.15 Errors on the quantity of interest p@,,,,. The red (resp. blue) color is used for zones where 
the quantity of interest is overestimated (resp. underestimated) [9] 
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Chapter 6 A) 
Applications and Extensions: A Survey cieca; 
of Literature 


This chapter contains a literature survey of the work published by the authors in the 
timeframe of their collaboration, where the concept presented in this book have been 
applied to real-life industrial settings, and new methodologies have been developed. 

The listed contributions are grouped in the following themes: linear manifold 
learning in Sect. 6.1, nonlinear dimensionality reduction via auto-encoder in Sect. 6.2, 
piecewise linear dimensionality reduction via dictionary-based ROM-nets in Sect. 6.3 
and manifold learning of physics problems assisted by black-box regressors in 
Sect. 6.4. 


6.1 Linear Manifold Learning 


A priori hyper-reduction method: an adaptive approach Model reduction meth- 
ods are usually based on offline preliminar simulations to build the shape functions 
of the reduced order model (ROM) before the computation of the reduced state vari- 
ables. They are a posteriori approaches. Most of the time these offline computations 
are as complex as the simulation which we want to simplify by the ROM. The a priori 
reduction method proposed in [17] avoids such preliminary computations. It is an a 
priori approach based on the analysis of some state evolutions, such that all the state 
evolutions needed to perform the model reduction are described by an approximate 
ROM. The ROM and the state evolution are simultaneously improved by the method, 
thanks to an adaptive strategy. An initial set of known shape functions can be used 
to define the ROM to adapt, but it is not necessary. The adaptive procedure includes 
extensions of the subspace spanned by the shape functions of the ROM and selec- 
tions of the most relevant shape functions in order to represent the state evolution. 
The hyper-reduction is achieved by selecting a part of the integration points of the 
finite element model to forecast the evolution of the reduced state variables. In this 
method, both the number of degrees of freedom and the number of integration points 
are reduced. 
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A nonintrusive distributed reduced-order modeling framework for nonlinear 
structural mechanics—Application to elastoviscoplastic computations In [8] is 
proposed a framework that constructs reduced-order models for nonlinear structural 
mechanics in a nonintrusive fashion and can handle large-scale simulations. Three 
steps are carried out: (i) the production of high-fidelity solutions by commercial soft- 
ware, (ii) the offline stage of the model reduction, and (iii) the online stage where the 
reduced-order model is exploited. The nonintrusivity assumes that only the displace- 
ment field solution is known, and the proposed framework carries out operations on 
these simulation data during the offline phase. The compatibility with a new com- 
mercial code only needs the implementation of a routine converting the discretized 
solution into the used data format. The nonintrusive capabilities of the framework 
are demonstrated on numerical experiments using commercial versions of Z-set and 
Ansys Mechanical. The nonlinear constitutive equations are evaluated by using an 
external plugin. The large-scale simulations are handled using domain decomposi- 
tion and parallel computing with distributed memory. The features and performances 
of the framework are evaluated on two numerical applications involving elastovis- 
coplastic materials: the second one involves a model of high-pressure blade, where 
the framework is used to extrapolate cyclic loadings in 6.5 h, whereas the reference 
high-fidelity computation would take 9.5 days. 


Fast computation of transient thermal profiling of high-pressure compressors 
under non-parametrized boundary conditions variability In [9], a transient ther- 
mal problem is considered, with a nonlinear term coming from the radiation boundary 
condition and a nonparametrized variability in the form complex scenarios for the 
initial condition and the convection coefficients and external temperatures. A poste- 
riori reduced order modeling by snapshot Proper Orthogonal Decomposition is used. 
To treat the nonlinearity, hyperreduction is required since precomputing the poly- 
nomial nonlinearities becomes too expensive for the radiation term. The Empirical 
Cubature Method, originally proposed for nonlinear structural mechanics, is derived 
for these equations. The method is applied to the design of high-pressure compres- 
sors for civilian aircraft engines, where a fast evaluation of the solution temperature 
is required when testing new configurations. When using in the reduced solver the 
same model as the one from the high-fidelity code, the approximation is very accu- 
rate. However, when using a commercial code to generate the high-fidelity data, 
where the implementation of the model and solver is unknown, the reduced model is 
less accurate but still within engineering tolerances. Hence, the regularizing property 
of reduced order models, together with a nonintrusive approach, enables the use of 
commercial software to generate the data, even under some degree of uncertainty in 
the proprietary model or solver of the commercial software. 


Time Stable Reduced Order Modeling by an Enhanced Reduced Order Basis 
of the Turbulent and Incompressible 3D Navier-Stokes Equations In [3], the 
problem of constructing a time stable reduced order model of the 3D turbulent and 
incompressible Navier-Stokes equations is considered. The lack of stability asso- 
ciated with the order reduction methods of the Navier-Stokes equations is a well- 
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known problem and, in general, it is very difficult to account for different scales of a 
turbulent flow in the same reduced space. To remedy this problem, a new stabiliza- 
tion technique is proposed, based on an a priori enrichment of the classical proper 
orthogonal decomposition (POD) modes with dissipative modes associated with the 
gradient of the velocity fields. The main idea is to be able to do an a priori analysis of 
different modes in order to arrange a POD basis in a different way, which is defined 
by the enforcement of the energetic dissipative modes within the first orders of the 
reduced order basis. This enables the modeling of the production and dissipation 
of the turbulent kinetic energy (TKE) in a separate fashion within the high ranked 
new velocity modes, hence to ensure good stability of the reduced order model. The 
importance of this a priori enrichment of the reduced basis is illustrated on a typical 
aeronautical injector with Reynolds number of 45,000. This order reduction tech- 
nique is able to recover large scale features for very long integration times (25 ms in 
the present case). Moreover, the reduced order model exhibits periodic fluctuations 
with a period of 2.2 ms corresponding to the time scale of the precessing vortex core 
(PVC) associated with this test case. 


An updated Gappy-POD to capture non-parameterized geometrical variation 
in fluid dynamics problems In [4], a method is proposed to fill the gap within an 
incomplete turbulent and incompressible data field, while satisfying the topological 
and intensity changes of the fluid flow after a non-parameterized geometrical varia- 
tion in the fluid domain. A single baseline large eddy simulation (LES) is assumed 
to be performed prior geometrical variations. The method enhances the Gappy-POD 
method proposed by Everson and Sirovich in 1995, in the case where the given set of 
empirical eigenfunctions is not sufficient and is not interpolant for the recovering of 
the modal coefficients for each Gappy snapshot by a least squares procedure. This is 
typically the case when we introduce non-parameterized geometrical modifications 
in the fluid domain. Here, after the baseline simulation, additional solutions of the 
incompressible Navier-Stokes equations are solely performed over a restricted fluid 
domain, that contains the geometrical modifications. These local LESs, called hybrid 
simulations, are performed by using the immersed boundary technique, which uses 
of a fluid boundary condition and the baseline velocity field. Then, the POD modes 
are updated using a local modification of the baseline POD modes in the restricted 
fluid domain. Furthermore, a physical correction of the latter enhanced Gappy-POD 
modal coefficients is obtained thanks to a Galerkin projection of the Navier-Stokes 
equations upon the new modes of the available data. This enhancement procedure 
on the global velocity reconstruction by the physical constraint was tested on a 3D 
semi-industrial test case of a typical aeronautical injection system and, a 2D laminar 
and unsteady incompressible test case. The speed-up relative to this new technique is 
equal to 100, which allows the fast exploration of two new designs of the aeronautical 
injection system. 
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6.2 Nonlinear Dimensionality Reduction via Auto-Encoder 


Data-Targeted Prior Distribution for Variational AutoEncoder In [1], Variational 
AutoEncoders (VAE) are used to study unsteady and compressible fluid flows in air- 
craft engines. Inferential methods enable the computation of sharp approximations of 
the posterior probability of the parameters of the VAE using the transient dynamics 
of the training velocity fields, and the generatation of plausible velocity fields. An 
important application is the initialization of such transient simulations. It is known 
by the Bayes theorem that the choice of the prior distribution is very important for 
the computation of the posterior probability, proportional to the product of likelihood 
with the prior probability. A new inference model is proposed, based on a new prior 
defined by the density estimate with the realizations of the kernel proper orthogonal 
decomposition coefficients of the available training data. It is illustrated that this 
inference model improves the results obtained with the usual standard normal prior 
distribution. This new generative approach can also be seen as an improvement of 
the kernel proper orthogonal decomposition method, for which we do not usually 
have a robust technique for expressing the pre-image in the input physical space of 
the stochastic reduced field in the feature high-dimensional space with a kernel inner 
product. 


A Bayesian Nonlinear Reduced Order Modeling Using Variational AutoEn- 
coders In [2], a new nonlinear projection based model reduction using convolutional 
VAEs. This framework is applied on transient incompressible flows. The accuracy 
is obtained thanks to the expression of the velocity and pressure fields in a nonlinear 
manifold maximising the likelihood on pre-computed data in the offline stage. A 
confidence interval is obtained for each time instant thanks to the definition of the 
reduced dynamic coefficients as independent random variables for which the pos- 
terior probability given the offline data is known. The parameters of the nonlinear 
manifold are optimized as the ones of the decoder layers of an autoencoder. The 
parameters of the conditional posterior probability of the reduced coefficients are the 
ones of the encoder layers of the same autoencoder. This reduced-order model is not 
a regression model over the offline pre-computed data. The numerical resolution of 
the ROM is based on the Chorin projection method. This new nonlinear projection- 
based reduced order modeling is applied to a 2D Karman Vortex street flow and a 
3D incompressible and unsteady flow in an aeronautical injection system. 


Deep multimodal autoencoder for crack criticality assessment In continuum 
mechanics, the prediction of defect harmfulness requires to solve approximately 
partial differential equations with given boundary conditions. In [14], boundary 
conditions are learnt for tight local volumes (TLV) surrounding cracks in three- 
dimensional volumes. A nonparametric data-driven approach is used to define the 
space of defects, by considering defects observed via X-Ray computed tomography. 
The dimension of the ambient space for the observed images of defects is huge. A 
nonlinear dimensionality reduction scheme is proposed in order to train a reduced 
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latent space for both the morphology of defects and their local mechanical effects in 
the TLV. A multimodal autoencoder enables to mix morphological and mechanical 
data. It contains a single latent space, termed mechanical latent space. But this latent 
space is fed by two encoders. One is related to the images of defects and the other to 
mechanical fields in the TLV. The latent variables are input variables for a geomet- 
rical decoder and for a mechanical decoder. In this work, mechanical variables are 
displacement fields. The autoencoder on mechanical variables enables projection- 
based model order reduction as proposed in the study of Lee and Carlberg. The main 
novelty of this work is a submodeling approach assisted by artificial intelligence. 
Here, for defect images in the test set, Dirichlet boundary conditions are applied to 
TLV. These boundary conditions are forecasted by the mechanical decoder with a 
latent vector predicted by the morphological encoder. For that purpose, a mapping 
is trained to convert morphological latent variables into mechanical latent variables, 
denoted “direct mapping”. An “inverse mapping” is also trained for error estima- 
tion with respect to morphological predictions. Errors on mechanical predictions are 
close to 5% with simulation speed-up ranging for 3 to 120. It is illustrated that latent 
variables forecasted by the images of defects are prone to a better understanding of 
the predictions. 


6.3 Piecewise Linear Dimensionality Reduction via 
Dictionary-Based ROM-Nets 


Model order reduction assisted by deep neural networks (ROM-net) In [12], a 
general framework is proposed for projection-based model order reduction assisted 
by deep neural networks. The proposed methodology, called ROM-net, consists in 
using deep learning techniques to adapt the reduced-order model to a stochastic 
input tensor whose nonparametrized variabilities strongly influence the quantities 
of interest for a given physics problem. In particular, the concept of dictionary- 
based ROM-nets is introduced, where deep neural networks recommend a suitable 
local reduced-order model from a dictionary. The dictionary of local reduced-order 
models is constructed from a clustering of simplified simulations enabling the iden- 
tification of the subspaces in which the solutions evolve for different input tensors. 
The training examples are represented by points on a Grassmann manifold, on which 
distances are computed for clustering. This methodology is applied to an anisother- 
mal elastoplastic problem in structural mechanics, where the damage field depends 
on a random temperature field. When using deep neural networks, the selection of 
the best reduced-order model for a given thermal loading is 60 times faster than when 
following the clustering procedure used in the training phase. 


Mechanical dissimilarity of defects in welded joints via Grassmann manifold 
and machine learning Assessing the harmfulness of defects based on images is 
becoming more and more common in industry. At present, these defects can be 
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inserted in digital twins that aim to replicate in a mechanical model what is observed 
on a component so that an image-based diagnosis can be further conducted. How- 
ever, the variety of defects, the complexity of their shape, and the computational 
complexity of finite element models related to their digital twin make this kind of 
diagnosis too slow for any practical application. In [18], a classification of observed 
defects enables the definition of a dictionary of digital twins. These digital twins 
prove to be representative of model-reduction purposes while preserving an accept- 
able accuracy for stress prediction. Nonsupervised machine learning is used for both 
the classification issue and the construction of reduced digital twins. The dictionary 
items are medoids found by a k-medoids clustering algorithm. Medoids are assumed 
to be well distributed in the training dataset according to a metric or a dissimilarity 
measurement. A new dissimilarity measurement between defects is proposed. It is 
theoretically founded according to approximation errors in hyper-reduced predic- 
tions. In doing so, defect classes are defined according to their mechanical effect 
and not directly according to their morphology. In practice, each defect in the train- 
ing dataset is encoded as a point on a Grassmann manifold. This methodology is 
evaluated through a test set of observed defects totally different from the training 
dataset of defects used to compute the dictionary of digital twins. The most appro- 
priate item in the dictionary for model reduction is selected according to an error 
indicator related to the hyper-reduced prediction of stresses. No plasticity effect is 
considered here (merely isotropic elastic materials), which is a strong assumption 
but which is not critical for the purpose of this work. In spite of the large variety of 
defects, accurate predictions of stresses for most of defects in the test set are obtained. 


Multimodal data augmentation for digital twining assisted by artificial intelli- 
gence in mechanics of materials Digital twins in the mechanics of materials usually 
involve multimodal data in the sense that an instance of a mechanical component has 
both experimental and simulated data. These simulations aim not only to replicate 
experimental observations but also to extend the data. Whether spatially, temporally, 
or functionally, augmentation is needed for various possible uses of the components to 
improve the predictions of mechanical behavior. Related multimodal data are scarce, 
high-dimensional and a physics-based causality relation exists between observational 
and simulated data. In [5], a data augmentation scheme coupled with data pruning is 
proposed, in order to limit memory requirements for high-dimensional augmented 
data. This augmentation is desirable for digital twining assisted by artificial intelli- 
gence when performing nonlinear model reduction. Here, data augmentation aims 
at preserving similarities in terms of the validity domain of reduced digital twins. 
A specimen subjected to a mechanical test at high temperature is considered, where 
the as-manufactured geometry may impact the lifetime of the component. Hence, an 
instance is represented by a digital twin that includes 3D X-Ray tomography data of 
the specimen, the related finite element mesh, and the finite element predictions of 
thermo-mechanical variables at several time steps. There is, thus, for each specimen, 
geometrical and mechanical information. Multimodal data, which couple different 
representation modalities together, are hard to collect, and annotating them requires 
a significant effort. Thus, the analysis of multimodal data generally suffers from the 
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problem of data scarcity. The proposed data augmentation scheme aims at training 
a recommending system that recognizes a category of data available in a training set 
that has already been fully analyzed by using high-fidelity models. Such a recom- 
mending system enables the use of a ROM-net for fast lifetime assessment via local 
reduced-order models. 


Real-Time Data Assimilation in Welding Operations Using Thermal Imaging 
and Accelerated Digital Twinning Welding operations may be subjected to dif- 
ferent types of defects when the process is not properly controlled and most defect 
detection is done a posteriori. The mechanical variables that are at the origin of these 
imperfections are often not observable in situ. In [16], an offline/online data assimila- 
tion approach is proposed, that allows for joint parameter and state estimations based 
on local probabilistic surrogate models and thermal imaging in real-time. Offline, 
the surrogate models are built from a high-fidelity thermomechanical finite element 
parametric study of the weld. The online estimations are obtained by conditioning 
the local models by the observed temperature and known operational parameters, 
thus fusing high-fidelity simulation data and experimental measurements. 


Mechanical assessment of defects in welded joints: morphological classifica- 
tion and data augmentation In [15], a methodology is developed for classifying 
defects based on their morphology and induced mechanical response. The proposed 
approach is fairly general and relies on morphological operators and spherical har- 
monic decomposition as a way to characterize the geometry of the pores, and on 
the Grassman distance evaluated on FFT-based computations, for the predicted elas- 
tic response. The approach is implemented and detailed on a set of trapped gas 
pores observed in X-ray tomography of welded joints, that significantly alter the 
mechanical reliability of these materials. The space of morphological and mechan- 
ical responses is first partitioned into clusters using the “k-medoids” criterion and 
associated distance functions. Second, multiple-layer perceptron neuronal networks 
are used to associate a defect and corresponding morphological representation to 
its mechanical response. It is found that the method provides accurate mechanical 
predictions if the training data contains a sufficient number of defects representing 
each mechanical class. To do so, the original set of defects is supplemented by data 
augmentation techniques. Artificially-generated pore shapes are obtained using the 
spherical harmonic decomposition and a singular value decomposition performed on 
the pores signed distance transform. 


Data Augmentation and Feature Selection for Automatic Model Recommenda- 
tion in Computational Physics Classification algorithms have recently found appli- 
cations in computational physics for the selection of numerical methods or models 
adapted to the environment and the state of the physical system. For such classifi- 
cation tasks, labeled training data come from numerical simulations and generally 
correspond to physical fields discretized on a mesh. Three challenging difficulties 
arise: the lack of training data, their high dimensionality, and the non-applicability 
of common data augmentation techniques to physics data. In [13], two algorithms 
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are introduced to address these issues: one for dimensionality reduction via feature 
selection, and one for data augmentation. These algorithms are combined with a wide 
variety of classifiers for their evaluation. When combined with a stacking ensem- 
ble made of six multilayer perceptrons and a ridge logistic regression, they enable 
reaching an accuracy of 90% on the considered classification problem of nonlinear 
structural mechanics. 


Physics-informed cluster analysis and a priori efficiency criterion for the con- 
struction of local reduced-order bases Nonlinear model order reduction has opened 
the door to parameter optimization and uncertainty quantification in complex physics 
problems governed by nonlinear equations. In particular, the computational cost 
of solving these equations can be reduced by means of local reduced-order bases. 
In [11], the benefits of a physics-informed cluster analysis for the construction of 
cluster-specific reduced-order bases are examined. The choice of the dissimilarity 
measure for clustering is fundamental and highly affects the performances of the 
local reduced-order bases. It is shown that clustering with an angle-based dissimilar- 
ity on simulation data efficiently decreases the intra-cluster Kolmogorov N-width. 
Additionally, an a priori efficiency criterion is introduced to assess the relevance of 
ROM-nets. This criterion also provides engineers with a very practical method for 
ROM-nets’ hyperparameters calibration under constrained computational costs for 
the training phase. On five different physics problems, the physics-informed cluster- 
ing strategy significantly outperforms classic strategies for the construction of local 
reduced-order bases in terms of projection errors. 


6.4 Extension: Manifold Learning of Physics Problems 
Assisted by Black-Box Regressors 


A priori compression of convolutional neural networks for wave simulators Con- 
volutional Neural Networks (CNNs) are seeing widespread use in a variety of fields, 
including image classification, facial and object recognition, medical imaging anal- 
ysis, and many more. In addition, there are applications such as physics-informed 
simulators in which accurate forecasts in real time with a minimal lag are required. 
The current neural network designs include millions of parameters, which makes it 
difficult to install such complex models on devices that have limited memory. Com- 
pression techniques might be able to resolve these issues by decreasing the size of 
CNN models that are created by reducing the number of parameters that contribute 
to the complexity of the models. In [7], a compressed tensor format of convolutional 
layer is proposed a priori, before the training of the neural network, for finite element 
(FE) predictions of physical data. 3-way kernels or 2-way kernels in convolutional 
layers are replaced by one-way fiters. The overfitting phenomena will be reduced 
also. The time needed to make predictions or time required for training using the 
original CNN model would be cut significantly if there were fewer parameters to 
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deal with. The priori compressed models is validated on physical data from a FE 
model solving a 2D wave equation. Tn the considered application, the proposed con- 
volutional compression technique achieves equivalent performance in the prediction 
error as Classical convolutional layers with fewer trainable parameters (around 20%) 
and lower memory footprint. 


Accelerated uncertainty quantification in impact simulations using genera- 
tive adversarial networks and submodeling The analysis of parametric and non- 
parametric uncertainties of very large dynamical systems requires the construction of 
a stochastic model of said system. Linear approaches relying on random matrix theory 
and principal component analysis can be used when systems undergo low-frequency 
vibrations. In the case of fast dynamics and wave propagation, a random generator 
of boundary conditions for fast submodels by using machine learning is investigated 
in [6]. It is illustrated that the use of non-linear techniques in machine learning and 
data-driven methods is highly relevant. Physics-Informed Neural Networks (PINNs) 
are a possible choice for a data-driven method to replace linear modal analysis. An 
architecture that supports a random component is necessary for the construction of 
the stochastic model of the physical system for non-parametric uncertainties, since 
the goal is to learn the underlying probabilistic distribution of uncertainty in the data. 
Generative Adversarial Networks (GANs) are suited for such applications, where the 
Wasserstein-GAN with gradient penalty variant offers improved convergence results 
for the considered problem. The objective of this approach is to train a GAN on data 
from a finite element method code so as to extract stochastic boundary conditions 
for faster finite element predictions on a submodel. The submodel and the training 
data have both the same geometrical support. It is a zone of interest for uncertainty 
quantification and relevant to engineering purposes. In the exploitation phase, the 
framework can be viewed as a randomized and parametrized simulation generator 
on the submodel, which can be used as a Monte Carlo estimator. 


MMGBP: a Mesh Morphing Gaussian Process-based machine learning method 
for regression of physical problems under non-parameterized geometrical vari- 
ability When learning simulations for modeling physical phenomena in industrial 
designs, geometrical variabilities are of prime interest. While classical regression 
techniques prove effective for parameterized geometries, practical scenarios often 
involve the absence of shape parametrization during the inference stage, providing 
only mesh discretizations as available data. Learning simulations from such mesh- 
based representations poses significant challenges, with recent advances relying 
heavily on deep graph neural networks to overcome the limitations of conventional 
machine learning approaches. Despite their promising results, graph neural networks 
exhibit certain drawbacks, including their dependency on extensive datasets and limi- 
tations in providing built-in predictive uncertainties or handling large meshes. In [10], 
a machine learning method that do not rely on graph neural networks is proposed. 
Complex geometrical shapes and variations with fixed topology are dealt with using 
well-known mesh morphing onto acommon support, combined with classical dimen- 
sionality reduction techniques and Gaussian processes. The proposed methodology 
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can easily deal with large meshes without the need for explicit shape parameteriza- 
tion and provides crucial predictive uncertainties, which are essential for informed 
decision-making. In the considered numerical experiments, the proposed method 
is competitive with respect to existing graph neural networks, regarding training 
efficiency and accuracy of the predictions. 
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